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ABSTRACT 

This  report  is  a  siirvey  of  the  function  approximation  routines 
available  for  the  BRLESC  computer  at  Aberdeen  Proving  Ground  as  of 
April  1969*  It  includes  a  description  of  each  routine  and  some  general 
discussion. 
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I.  INTRODUCTION 


The  approximation  of  a  function  by  another  function  is  a  frequent 

computing  problem.  The  primary  purpose  of  this  report  is  to  list  the 

1* 

readily  available  routines  of  this  type  for  the  BRIESC  at  Aberdeen 
Proving  Ground^  as  of  April  I969.  The  majority  of  these  routines,  and 
descriptions  of  them,  are  available  through  the  Computer  Support 
Division,  Systems  Programming  (Building  528,  Room  215).  Anyone  who  is 
considering  Curve  fitting  is  strongly  urged  to  examine  the  existing 
programs  before  attempting  to  code  his  own. 

In  nearly  all  of  these  routines  the  original  function  is 
represented  by  a  table  of  N  data  points,  (X^,y^),  i=l,2, ...,N,  where  y^ 
is  the  dependent  variable  (the  function  value  we  wish  to  approximate) 

and  X.  is  either  a  single  dependent  variable,  x.,  or  a  vector  of 

1  T  ^ 

variables,  X.  =  (x. ^ ,x. . ,X.„)  .  The  immediate  goal  of  most  of  the 
^  1  '  il’  i2^  ’  iN' 

routines  is  to  find  the  "best"  values  for  M  parameters  A  = 

in  an  equation  F(X,A),  the  approximating  function,  which  has  been 

chosen  to  represent  the  table  of  data. 

In  all  these  routines,  the  parameters.  A,  are  chosen  so  as  to 
satisfy  some  prescribed  constraints  or  to  optimize  some  measure  of  the 
goodness  of  the  approximation.  Some  of  the  routines  do  both.  The  most 
popular  measure  is  the  variance  of  residuals  (the  sum  of  squares  of  the 
residuals)  over  the  set  of  data  being  fitted.  Fitting  with  this 
measure  is  called  the  "method  of  least  squares".  Most  numerical 
analysis  textbooks  have  some  discussion  of  the  principle  of  least 
squares  (e.g.  Hildebrand^).  Most  of  the  fitting  programs  for  BRLESC 
are  based  on  this  principle  which  states  that  the  "best"  approximation 
is  that  which  minimizes  the  sum  of  squares  of  the  residuals. 

The  only  approximation  routines  we  consider  are  for  interpolation 
or  curve  fitting.  We  will  not  discuss  the  approximations  used  for  any 
of  the  standard  functions  (SIN,  EXP,  etc.),  nor  will  we  consider  any 


References  are  listed  on  pages  UO  and  Ul. 
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integration  routines.  Undoubtedly,  some  useful  programs  were  overlooked, 

other  good  BRLESC  approximation  routines  may  be  completed  in  the  near 

p  5 

future,  and  FORTRA.N  versions  of  the  most  useful  FORAST  routines  will 
be  made.  (The  MULTIPLE  REGRESSION  and  NONLINEAR  LEAST  SQUARES  are  now 
available  in  FORTRAN.) 

In  this  report  we  first  list  the  FORTRAN  and  FORAST  function 
approximation  routines  which  were  located  and  give  a  brief  description 
of  them.  We  then  have  some  general  discussion  of  these  programs. 

Finally,  we  discuss  each  routine  in  a  little  more  detail,  give 
references,  point  out  some  difficulties,  and  indicate  some  advantages 
and  disadvantages  of  some  of  the  routines. 

Instructions  for  using  many  of  the  routines  are  included  in  the 
appendix . 


II.  LIST  OF  ROUTINES 

A.  Standard  FORTRAN  and  ( FORAST)  Subroutines 
Interpolation  Subroutine 
DVDINT  (D.D.IN) 


Least  Squaxes  Subroutines 


MATINV,  FNEQS  (S.N.E..  F.N.E.,  etc.) 
GENLSQ  (G.L.SQ) 

POLYLS  (P.L.SQ) 


Function  Minimizing  Subroutines 

FNMIN  (FN.MIN) 

FDMIN  (FD.MIN) 


B.  FORAST  Fitting  Programs 

LEAST  SQUARES  PROGRAM 

MULTIPLE  REGRESSION  (FORTRAN  program  available) 
NONLINEAR  LEAST  SQUARES  (FORTRAN  program  available) 
PIECEWISE  QUARTIC  FTT 
CUBIC  SPLINE 

LEAST  SQUARES  CUBIC  SPLINE 
POLYGONAL  CURVE 
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C.  Other  FORTRAN  Routines 

CONSTRAINED  NONLINEAR  LEAST  SQUARES  (Program) 

NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA 

NLPROG  (Subroutine  to  minimize  constrained  functions) 

The  "standard"  FORTRAN  subroutines  in  the  above  list  are  available 
on  cards  from  Systems  Programming,  Computer  Support  Division.  Each  of 
these  subroutines  is  equivalent  to  some  FORAST  Subroutine  on  the  FORAST 
compiler  tape.  Three  of  the  FORAST  programs,  LEAST  SQUARES  PROGRAM, 
MULTIPLE  REGRESSION,  and  NONLINEAR  LEAST  SQUARES,  are  also  available 
from  Systems  Programming.  Using  routines  from  this  source  has  several 
advantages:  The  routines  are  available  in  a  standard  usable  form. 
Descriptions  aj^e  available.  The  routines  are  completely  checked  out. 

And  advice  and  assistance  are  available  if  needed.  The  other  FORAST 
programs  and  the  "other"  FORTRAN  routines  have  special  features  which 
may  compensate  for  their  defects  in  availability,  standardization,  etc. 

We  have  classified  these  routines  as  interpolating  routines, 
fitting  routines,  least  squares  routines,  and  function  minimizing 
routines.  Such  a  division  must  be  arbitrary.  We  have  elected  to  call 
the  method  "interpolation"  if  the  principal  result  is  the  value  of  the 
approximating  function  at  a  prescribed  value  of  the  independent  variable. 
If  the  main  output  of  the  routine  is  an  equation  or  the  values  of 
parameters  in  a  prescribed  equation,  we  call  it  a  "fitting"  routine. 

Those  fitting  procedures  which  use  the  principle  of  least  squares  are 
referred  to  as  "least  squares"  routines. 

It  seems  to  be  general  practice  to  drop  the  adjective  "linear"  in 
discussion  of  linear  least  squares  and  to  use  the  adjective  "nonlinear" 
when  referring  to  least  squares  procedures  which  are  not  necessarily 
linear  in  the  unknown  parameters.  This  practice  has  been  followed 
here  in  naming  and  discussing  the  routines.  To  further  complicate  these 
titles,  the  adjectives  "general"  and  "polynomial"  are  used  to  denote 
the  type  of  terms  permitted  in  linear  least  squares. 
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The  function  minimizing  routines  are  separated  from  the  other 
routines  because  the  user  is  free  to  choose  the  measure  of  the  goodness 
of  the  approximation.  The  user  pays  for  this  freedom  with  coding  and 
some  loss  of  confidence  in  the  results. 


DVDINT  (D.D.IM  in  FORAST) 

This  standard  subroutine  uses  divided  differences  to  approximate 
the  value  of  y  corresponding  to  x  from  a  table  (x^,y^),  i=l,2, ...,N, 
ordered  on  x.  The  approximating  function  is  an  M-th  degree  polynomial 
through  M  +  1  points  in  the  neighborhood  of  x. 

MATINV,  FNEQS  (S.N.E.  ,  F.N.E.  ,  etc,  in  FQFtAST) 

The  MATINV  subroutine  is  the  standard  FORTRAN  subroutine  for 
solving  a  system  of  linear  equations.  The  FNEQS  subroutine  may  be  used 
to  help  convert  raw  data  to  the  normal  equations  of  least  squares  and 
the  MATINV  subroutine  may  be  used  to  solve  this  system.  In  FORAST,  the 
F.N.E.  may  be  used  to  help  form  the  normal  equations,  in  the  form  of  an 
augmented  upper  triangular  matrix,  and  the  SY.SNE  or  SY.INV  used  to 
solve  the  system.  One  can  also  use  the  F.O.MAT  to  form  the  full 
rectangular  augmented  matrix  from  the  triangular  form  and  use  the  S.N.E. 
or  MAT.INV.  to  solve  the  system. 


It  is  generally  better  to  use  one  of  the  next  two  subroutines  for 
linear  least  squares. 

GENLSQ  (G.L.SQ  in  FORAST) 

Starting  from  a  table  of  data  (y^,cp^(X^)  ,cp2(X^) , . .  .  ,cpj^(X^) ) , 
i=l,2,...,N,  this  routine  carries  out  a  linear  least  squares  fit.  It 


IS  a 


finds  A  =  (a^,a2,...,a^)  so  that  U(A)  = 

minimum  with  respect  to  the  a.'s.  This  routine  will  also  compute  some 
information  to  help  the  user  analyze  the  fit:  the  root -mean-square 


error 


1 

,  (U(A)/(N-M)  )^;  the  residuals,  R.  =  y.  -)  a.cp.(X.);  a., 

1  1  ^J-l  J  J  1  J 


the 


estimates  of  the  error  in  a.:  and  t.  =  a. /a.,  the  estimates  of  the 

J  J  J  J 

significance  of  a . . 

J 
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The  FOPtAST  subroutine  allows  more  flexibility  in  the  storage  of  the 
input  data  and  permits  more  flexibility  in  weighting  or  eliminating  data 
points . 

POLYLS  (P.L.SQ  in  FORAST) 

This  is  simply  a  special  case  of  GENLSQ  (G.L.SQ)  for  polynomial 
least  squares  with  one  independent  variable.  The  input  is  N  pairs  of 
numbers  (x^,y^).  The  program  minimizes  the  equation 

U(A)  =  y^_-^  (y^.  -  The  FORAST  program  has  an  option  for 

dropping  some  of  the  powers  of  x. 

FNMIU  (FN.MIN  in  FORAST) 

This  is  a  standard  subroutine  for  finding  the  minimum  of  a  function 
of  several  variables  without  using  derivatives.  The  least  squares 

2  ■g' 

programs  we  have  Just  discussed  find  the  minimum  of  U(A)  =  (  Z.i=i^)  ' 

where  the  are  linear  in  a^^a^^ • • • This  minimizing  program  could 

be  used  to  minimize  U(A)  if  the  R.  are  not  linear  in  the  a.'s  or  if 

1  J 

some  other  measure  of  the  goodness  of  the  fit  were  used.  The  user  must 
supply  programming  to  evaluate  U(A). 

FDMIN  (FD.MIN  in  FORAST) 

This  is  another  function  minimizing  routine.  It  uses  the  first 

partial  derivatives  of  the  function  to  be  minimized.  The  user  must 

supply  programming  which  evaluates  U(A)  and  9U(A)/9a.,  j=l,2,...jM. 

J 

LEAST  SQUARES  PROGRAM 

This  complete  FORAST  program  reads  control  cards  that  specify  the 
type  and  form  of  a  linear  least  squares  fit,  reads  the  data,  computes 
the  fit,  and  prints  the  coefficients,  the  approximate  function  value  at 
each  point,  the  corresponding  residuals,  the  root-mean-square  error,  and 
the  "sigma  and  t"  error  indicators. 
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MULTIPLE  REGRESSION  PRCX}RAM 

This  complete  FORAST  program  is  similar  to  the  previous  program  in 
set  up  and  control.  The  very  important  difference  is  that  the  given 
approximating  equation  is  a  listing  of  possible  terms  to  be  included  in 
a  linear  least  squares  fit.  This  routine  selects  an  equation  that  has 
only  significant  terms  and  then  carries  out  a  least  squares  fit  to  find 
the  best  coefficients  for  this  equation.  In  addition  to  output  like  the 
LEAST  SQUARES  PROGRAM,  this  program  prints  a  running  account  of  the 
terms  added  to  or  dropped  from  the  tentative  equation  and  prints  the 
terms  in  the  final  equation. 

NONLINEAR  LEAST  SQUARES 

This  complete  FORAST  program  is  a  program  for  finding  the  "best" 
parameters  in  the  least  squares  sense  for  an  approximating  function 
which  is  not  linear  in  all  the  parameters.  If  the  function  F(X,A)  is 
used  to  approximate  y,  the  best  approximation  in  the  least  squares 

sense  occurs  when  U(A)  =  "  F(X^,A))^  is  a  minimimi.  F(X,A)  is 

approximated  with  the  linear  terms  of  Taylor’s  series  (i.e., 

—  — 

F(X^,A)  =  F(X^,A)  +  2_^j_^Aa^dF^/9aj ,  where  A  is  a  good  initial 

approximation  to  A,  and  dF./ba.  is  the  partial  derivative  of  F(X,A) 

evaluated  for  A  =  A  and  X  =  X^).  The  standard  procedure  for  linear 

least  squares  is  used  to  find  Aa.,  j=l,2, ...,M,  so  that 

J 

U(AA)  =).  ,(y,  -  F(X,A)  -  )  _Aa,dF7da,)^  is  minimized.  Then  a., 

i-il-X  1  XjJ-X  J  1  j  j ' 

j=l,2,...,M,  is  replaced  by  a.  +  pAa.,  0  <  p  ^  1.  With  luck,  and  a 

J  J 

good  initial  guess,  this  method  will  converge  to  the  "best"  answer. 

This  program  requires  more  input  information  than  the  linear  LEAST 
SQUARES  PROGRAM,  In  addition  to  output  like  the  LEAST  SQUARES  PROGRAM, 
this  program  prints  information  about  each  iteration. 
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PIECEWISE  QUARTIC  FIT 


This  complete  FORAST  program  approximates  a  bivariate  table  of  data 
(Xi^Yi),  i=l,2,...,N,  with  a  series  of  connected  quartic  polynomials. 

As  many  points  as  possible  (without  exceeding  a  specified  least  squares 
error)  are  placed  in  each  succeeding  polynomial.  The  resulting 
approximating  fimction  is  continuous  and  has  continuous  first  derivatives. 

CUBIC  SPLINE 

This  complete  FORAST  program  fits  N  data  points  (x^,y^),  i=l,2,...,W, 
with  N  -  1  cubic  equations.  The  approximating  function  passes  through 
each  data  point,  is  continuous,  and  has  continuous  first  and  second 
derivatives . 

LEAST  SQUARES  CUBIC  SPLINE 

This  complete  FORAST  program  fits  N  data  points  (x^,y^),  t=l,2,...,N, 
with  a  prescribed  nijmber,  say  K,  cubic  equations.  The  user  specifies  the 
abscissas  for  K  +  1  "break  points".  The  program  finds  the  K  +  1 
corresponding  ordinates  so  that  the  cubic  spline  function  F(x),  through 

these  break  points,  makes  ^^(y^  -  F(x^))^  a  minimm. 

POLIGOWAL  CURVE 

It  is  sometimes  useful  to  approximate  bivariate  data  by  a  series  of 
connected  straight  lines.  This  program  uses  dynamic  programming  to 
locate  the  best  values  of  the  independent  variable  for  the  ends  of  a 
specified  number  of  straight  lines.  For  these  "break  points",  the 
approximating  polygonal  curve  is  the  best  in  the  least  squares  sense . 

CONSTRAINED  NONLINEAR  LEAST  SQUARES 

Sometimes  we  want  to  approximate  a  set  of  data  (X^,y^),  i=l,2,...,N, 
with  a  function  F(X,A)  and  at  the  same  time  satisfy  a  number,  say  K, 
inequality  (and/or  equality)  constraints  G,  (A)  S  0  (or  G  (A)  =  0), 

iC  iC 

k=l,2,...,K.  This  FORTRAN  program  finds  an  approximate  answer  to  this 
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problem  by  minimizing  the  function  U(A)  =  ^  ’  where 

f^i^)  =  -  F(X.,A)  for  i=l,2,...,N,  and  =  Gj^(A)  for 

unsatisfied  constraints  (fj^_^j^(A)  =  0  if  Gj^(A)  is  an  inequality 
constraint  greater  than  O). 

The  user  must  supply  three  subroutines:  FCODE  to  evaluate  f^(A), 
PCODE  to  evaluate  j=l,2,...,M,  and  SUEZ  to  set  up  some  initial 

conditions.  The  program  uses  a  combination  of  the  differential 
correction  method  (used  by  the  FORAST  NONLINEAR  LEAST  SQUARES  program) 
and  steepest  descent  to  minimize  U(A).  This  method  should  converge  in 
most  cases. 

NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA 

This  FORTRAN  routine  is  somewhat  different  than  the  other  least 

squares  routines.  The  data  points  consist  of  > 

m. ,  the  standard  error  of  weight  one;  and  R.,  the  correlation  matrix  of 
1  J 

cofactors.  The  object  of  the  program  is  to  find  "  ^i  "  ^i^ 

i=l,2,...,N,  and  A  =  so  thatU(|)  =  is 

minimized  and  F(X^  +  ^  i=l,2,...,N.  The  user  must  supply 

a  main  program  which  sets  up  the  input  and  controls  the  iterations,  and 
a  subroutine  to  compute  F(X^  +  5^jA),  and  the  partial  derivatives  of 
this  F  with  respect  to  k=l,2,...,K,  and  a.,  j=l,2,...,M. 

lA  j 

NLPROG 

This  is  a  function  minimizing  FORTRAN  subroutine  which  permits 
constraints.  If  U(A)  is  the  function  to  be  minimized  subject  to  the 
inequality  constraints  ^  0^  k=l,2,...,J,  and  the  equality 

constraints  G,  (A)  =  0,  k=J+l, J+2, . . . , J+K,  this  subroutine  is  designed  to 

2  ^ 

over  the  region 

Q  =  fA:Gj^(A)  s  0,  k=l,2,...,j]  for  p=p^,P2,  •  .  . ,  Pj^  with  p^  >  >  0 


minimize  P(A,p)  =  U(A)  +  p  ^^2 


If  the  original  estimate  for  A  is  not  in  Q, 


the  routine  will  first  try  to  find  an  initial  A  which  is  in  Q.  Actually, 

NLPROG  is  one  of  four  different  routines  depending  on  the  method  used 

to  minimize  P(A,p).  Two  of  these  are  direct  search  methods  which  do  not 

require  derivatives.  One  is  a  quasi-  Newton  method  which  uses  first 

derivatives  (the  partials  of  F  and  G  with  respect  to  a.).  The  other 

K  J 

routine  uses  Newton' s  method,  which  requires  first  and  second  partial 
derivatives . 


The  user  must  code  a  main  program  to  supply  the  initial  data  and 
control  parameters  for  NLPROG.  He  must  also  supply  a  subroutine  to 
evaluate  P  and  the  G^  and,  if  needed,  a  subroutine  to  evaluate  their 
derivatives . 


III.  GENERAL  DISCUSSION 
A.  Classification  of  Routines 

In  the  previous  section  we  characterized  the  BRLESC  function 
approximation  routines  as  FORAST  or  FORTRAN  routines,  as  standard  sub¬ 
routines  or  something  else,  and  as  "interpolating",  "fitting",  or 
"function  minimizing"  routines.  The  separation  between  FORAST  and 
FORTRAN  routines  is  necessary  and  definite.  The  distinction  between  the 
standard  subroutines  and  other  routines  is  also  useful  and  quite  clear. 
(The  standard  subroutines  are  all  on  the  FORAST  compiler  tape  and 
equivalent  FORTRAN  subroutines  are  available  on  cards.)  The  third 
classification  is  very  arbitreury.  The  broadest  definition  of  any  of 
the  three  terms  includes  all  the  routines.  An  "interpolating"  routine, 
by  our  definition,  is  one  that  produces  a  numerical  estimate  of  the 
dependent  variable  for  a  given  set  of  the  Independent  variables.  Thus, 
we  had  only  one  Interpolating  routine,  the  standard  divided  difference 
interpolation  routine.  It  is  convenient  to  classify  the  three 
"function  minimizing"  routines  separately  because  the  user  is  free  to 
choose  the  measure  of  "goodness"  of  the  fit.  All  the  remaining  "fitting" 
routines  except  the  CUBIC  SPLINE  use  the  principle  of  least  squares. 
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There  are  several  other  classifications  which  are  convenient  for 
some  purposes.  We  can  separate  the  bivariate  routines  and  the  routines 
that  can  be  used  with  more  than  two  variables.  We  can  also  classify 
the  routines  by  continuity,  smoothness,  or  by  the  measure  of  goodness 
of  the  fit.  Some  of  the  routines  use  or  permit  constraints.  These 
topics  will  be  taken  up  in  the  next  paragraphs. 

B.  Goodness  of  Approximation  and  Constraints 

Approximating  functions  generally  satisfy  some  set  of  constraints, 
optimize  some  measure  of  goodness  of  the  fit,  or  do  both.  The 
approximating  functions  for  two  of  our  BRLESC  routines,  the  DVDIWT  and 
the  CUBIC  SPLINE,  are  completely  determined  by  constraints.  For  the 
three  function  minimizing  routines,  the  function  to  be  minimized  is 
selected  by  the  user  (and  computed  by  a  user  coded  program).  All  the 
other  routines  use  the  principle  of  least  squares  and  minimize  the  sum 

2 

of  squares  of  the  residuals,  > 

t— il— X  1 


The  user  of  the  function  minimizing  routines  may  choose  the  sum  of 
squares  of  residuals,  U(A)  function  to  minimize.  (This 

is  equivalent  to  minimizing  so-called  "Euclidean"  norm.) 

Another  useful  choice  is  U(A)  =  Maximum|R^|,  i=l,2, ...,N,  the 
"Tchebycheff "  or  "uniform"  norm  used  to  find  the  minimax  or  Tchebycheff 

solution.  The  norm  U(A)  “2ji-ll^i^  another  frequently  mentioned 
choice.  The  form  of  the  approximating  function  is  selected  by  the  user. 


In  some  of  the  least  squares  routines  the  form  of  the  approximating 
function  is  specified  by  the  routine.  In  others  the  user  may  select 
this  function.  The  residuals  for  all  the  least  squares  routines,  except 
the  NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA  program,  are  the 
difference  between  the  given  value  of  the  dependent  variable  and  the 
approximated  value  (i.e.,  R^  =  y^-F(X^,A)).  The  measure  of  goodness  of 
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the  fit  is  the  unbiased  estimate  of  the  standard  deviation, 

ERMS  =  (2^^_^R^/(N-M)  )  ,  where  N  is  the  number  of  data  points  and  M  is 

the  number  of  coefficients  determined  by  the  routine . 

The  NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA  program  is  somewhat 
different.  The  N  data  points  consist  of  X.  =  (x. ,x.  ,  the 

measured  values  of  the  variables;  m^,  the  standard  error  of  weight  one; 
and  R. ,  the  correlation  matrix  of  cofactors.  The  coefficients 
A  =  (a^^a^,  •  •  •  ,a^)  and  the  correction  vectors  t=l,2,...,N,  are  found 

T  -1 

so  that  F(X^+5^,A)  =  0  for  all  i  and  U(|)  =  is  ^ 

minimum.  (u(^)/(N-M) is  used  as  the  measure  of  goodness  of  the  fit. 

The  FORTRAN  routine,  CONSTRAINED  NONLINEAR  LEAST  SQUARES  permits 
constraints.  For  this  program  the  user  must  compute  the  constraints 
Gj,  j=l,2, ...,J  as  well  as  the  residuals  R^,  i=l,2,...,N  and  their 
first  derivatives  so  he  may  define  them  as  he  wishes.  The  routine 

2  2  2  i 

minimizes  U*(A)  =  2ji-l^i  routine  prints  (2^^_^R^/(N-M) 

and  several  other  quantities  to  help  determine  the  validity  and  goodness 
of  the  approximation. 

C .  Bivariate  Routines 

Some  of  the  routines  are  designed  to  handle  many  variables  but  a 
few  are  restricted  to  two  variables,  a  dependent  variable  y,  and  an 
independent  variable  x.  All  of  these  approximate  y  by  a  polynomial  or 
set  of  polynomials  in  x.  The  routines  restricted  to  two  variables  are 
DVDINT  (D.D.IN),  POLYLS  (P.L.SQ),  PIECEWISE  QUARTIC  FIT,  CUBIC  SPLINE, 
LEAST  SQUARES  CUBIC  SPLINE,  and  POLYGONAL  CURVE.  The  first  of  these 
routines  is  polynomial  interpolation.  The  others  are  substitutes  for 
polynomial  interpolation  which  may  remove  or  reduce  the  main  objections 
to  interpolation. 

1.  Interpolation.  We  have  defined  an  interpolation  routine  as  one 
which  supplies  a  value  of  the  Independent  variable.  This  is  a  very 
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narrow  definition  but  it  describes  the  purpose  of  the  DVDINT  subroutine 
quite  well  From  experience  we  think  of  interpolation  as  a  process  by 
which  we  can  approximate  the  value  of  a  dependent  variable  between 
values  given  in  a  table .  If  we  include  approximation  for  the  given 
points  as  well,  we  are  describing  a  process  that  approximates  the 
tabular  function  by  some  other  function.  In  this  broader  sense,  all  the 
routines  discussed  in  this  report  are  interpolation  routines. 

Any  text  on  numerical  methods  devotes  some  space  to  interpolation. 
There  is  usually  some  mention  of  interpolation  involving  more  than  one 
independent  variable  and  some  discussion  about  using  various  functions, 
including  trigonometric  functions  to  interpolate  with;  but  the  major 
portion  of  the  discussion  is  restricted  to  interpolation  by  polynomials. 
Most  of  the  familiar  methods  of  interpolation  (e.g.,  Aitken’s  iterative 
interpolation,  Leigrange ' s  interpolation  formula,  Newton's  forward  and 
backward  difference  methods,  divided  difference,  etc.)  are  just  different 
ways  of  evaluating  an  M-th  degree  polynomial  through  M  +  1  points  near 
X  in  a  table  (x^,y^).  The  Newton's  interpolation  formula  with  divided 
differences  was  chosen  for  the  standard  interpolation  routine  (DVDINT 
in  FORTRAN,  D.D.IN  in  FORAST).  This  routine  is  easy  to  use,  efficient, 
and  accurate.  Unfortunately,  it  is  so  familiar  and  well  documented  that 
we  do  not  give  its  use  enough  thought,  let  alone  check  the  results. 
Polynomial  interpolation  has  some  drawbacks:  If  we  consider  the 
approximation  over  the  entire  range  of  the  table,  the  function  changes 
at  most  of  the  data  points  and  the  derivatives  fail  to  exist  at  these 
points,  even  though  the  function  itself  is  continuous.  A  high  degree 
polynomial  may  oscillate  wildly  and  errors  in  the  data  tend  to  augment 
this  oscillation.  This  is  pointed  out  in  the  Example  .  Finally,  it  is 
too  easy  to  overlook  gross  errors  in  the  data. 

2.  Other  Bivariate  Routines.  The  other  bivariate  routines  are,  in 
a  real  sense,  substitutes  for  polynomial  interpolation  which  remove  or 
reduce  one  or  more  of  the  objectionable  features  of  polynomial 
interpolation.  These  routines,  except  CUBIC  SPLINE,  use  least  squares 
and  consequently  reduce  the  effect  of  small  errors.  Some  of  these 
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routines  produce  all  the  residuals,  hence,  large  errors  in  the  data 
points  can  be  discovered.  The  other  routines  use  two  or  more  polynomials 
connected  at  points  we  call  "break  points". 

The  approximate  function  from  the  CUBIC  SPLINE  routine  passes 
through  each  data  point.  A  different  cubic  polynomial  is  used  between 
each  pair  of  data  points.  This  approximation  has  continuous  first  and 
second  derivatives. 

The  X  coordinates  of  the  break  points  for  the  LEAST  SQUARES  CUBIC 
SPLINE  routine  must  be  chosen  by  the  user.  The  routine  chooses  the 
corresponding  y  coordinates  so  that  the  cubic  spline  function  through 
the  break  points  produces  the  minimum  sum  of  squares  of  residuals.  This 
approximating  function  also  has  continuous  first  and  second  derivatives. 

The  break  points  for  the  POLYGONAL  CURVE  are  selected  by  the  program 
(The  user  must  state  the  number  of  them  to  be  used.)  There  are  no 
derivatives  at  these  break  points  since  the  straight  lines  which  meet 
there  have  different  slopes. 

The  PIECEWISE  QUARTIC  FIT  routine  selects  its  own  break  points 
from  the  data  points.  As  many  data  points  as  possible,  without  causing 
too  large  an  ERMS  error,  are  included  between  each  break  point.  The 
approximating  function  passes  through  each  of  these  break  points  and 
has  a  continuous  first  derivative  there. 

With  the  exception  of  the  divided  difference  Interpolation 
subroutine,  the  polynomial  least  squares  subroutine,  POLYLS,  is  the 
easiest  of  the  bivariate  routines  to  use.  If  a  good  fit  can  be  made 
with  a  low  degree  polynomial,  this  is  an  excellent  choice.  Unfortunately 
even  with  least  squares,  approximations  with  high  degree  polynomials  tend 
to  have  undesirable  oscillations. 

The  bivariate  routines,  with  the  exception  of  POLYLS  must  have  the 
data  ordered  on  the  independent  variable.  The  DVDINT  subroutine 
permits  either  increasing  or  decreasing  ordering.  The  other  routines 
assume  strictly  increasing  x^'s. 
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5-  Smoothness .  Although  it  is  possible  to  use  multivariate 

approximating  functions  which  have  discontinuities  or  other  disagreeable 

features,  the  usual  function  chosen  is  continuous  and  has  continuous 

derivatives  throughout  the  entire  range  of  definition.  This  is  not  as 

generally  true  for  bivariate  approximating  functions.  Approximation  by 

polynomial  interpolation  produces  an  approximating  function  which  does 

not  have  derivatives  at  most  of  the  data  points.  A  least  squares  fit  over 

the  entire  range  of  the  data  produces  an  approximating  function  with 

continuous  derivatives  of  all  orders;  but  fitting  with  a  low  degree 

polynomial  may  give  a  poor  approximation  at  the  data  points,  and  fitting 

with  a  high  degree  polynomial  may  give  good  results  at  the  data  points 

but  unbelievable  intermediate  values.  The  PIECEWISE  QUARTIC  FIT  program 

produces  an  approximation  that  has  continuous  first  derivatives  and 

should  not  oscillate  too  wildly.  The  CUBIC  SPLINE  and  the  LEAST  SQUARES 

CUBIC  SPLINE  programs  produce  an  approximating  function  that  has 

continuous  first  and  second  derivatives.  Spline  functions  have  received 

a  great  deal  of  attention  recently  because  of  their  "smoothness" 

properties.  In  particular,  the  cubic  spline  function,  say  S(x),  over 

2  2 

N  ^  5  data  points,  with  d  S/dx  =  0  at  x^  and  satisfies 


function  with  continuous  first  and  second  derivatives  which  satisfies 
f(x^)  =  y^;  i=l,2,...,N. 

4.  Checking.  A  graph  of  the  approximating  function  with  the 
original  data  points  superposed  gives  an  excellent  qualitative  check  of 
approximations  in  two  variables.  Gross  errors  in  the  input  data  are 
evident;  and  unusual  features  of  the  results,  such  as  excessive 
oscillation  or  poor  approximation  in  some  particular  region,  are  obvious. 
This  check  requires  extra  work,  but  it  is  worthwhile. 


2  2  2  2  2  2 

(d  S(x)/dx  )  dx  s  1  (d  f(x)/dx  )  dx  where  f(x)  is  any  continuous 
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D.  Least  Squares 


Most  of  the  fitting  programs  for  the  BRLESC  are  based  on  the 
principle  of  least  squares  applied  to  a  discrete  set  of  points 
i=l,2,...,N.  It  is  assiomed  there  is  some  function  Y(X)  for  which 
Y(X.)  =  y. .  (X.  may  be  a  single  variable  or  a  vector, 

X^  =  •)  ^he  function  Y{x)  is  to  be  approximated  by  a 

function  F(X,A)  where  A  represents  the  M  unknown  parameters  a^,a2, .••,a^ 
The  "best"  fit  in  the  least  squares  sense  (for  a  particular  fiinction 
F(X,A)  relative  to  a  non-negative  weighting  function  w(X),  over  the  set 
of  K  data  points)  is  achieved  when  A  is  found  such  that 


U(A)  =  )^.^^w(X.)(Y(X.) 


f(X.,A))^  is  a  minimum. 


For  this  to  be 


meaningful,  there  must  be  at  least  as  many  data  points  as  unknown 
parameters  (i.e.,  N  ^  M) .  We  will  assume  w(x)  =  1  in  what  follows. 
This  is  the  most  frequent  choice. 


1.  Linear  Least  Squares  Routines.  If  F(X,A),  the  approximating 
function,  is  linear  with  respect  to  the  M  unknown  parameters 

2 

a^,a2,...,aj^  we  can  write  U(A)  =  ^  ^ 

9.(X),  J=1,2,...,M,  are  M  suitable  functions  of  X.  (As  far  as  the  fit 
J 

is  concerned,  the  9.(X)  do  not  have  to  be  defined  except  at  the  N  data 

J 

points.  However,  if  the  result  is  to  be  useful  they  should  be  well 

behaved  and  fairly  easy  to  evaluate  for  all  X  in  the  range  of  the  fit.) 

The  function  U(A)  will  be  a  minimum  if  BU/Sa^^  =  0,  k=l,2,...,M.  This 

is  a  system  of  M  linear  equations,  called  "normal  equations",  in  the  M 

unknown  a . : 

J 


Zjj=i^j  ^=i\ 


k=l,2, . . . ,M. 


This  can  be  rewritten  in  matrix  notation  as  WA  =  V,  where  A  is  the  vector 
(a^a^, .  .  .  ,aj^)^,  W  is  the  symmetric  MxM  matrix  with  .l9.(X.)cPk(Xi) 
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T 

as  the  element  in  row  j  colijmn  k,  and  V  is  the  vector  •  ,Vj^) 

with  V,  =  )  .  ,cp,  (X.  )y..  (If  we  consider  cp.(X.  )  to  be  the  j-th  element 
k  Zji=l^k^  1  •'i  i 

in  the  i-th  row  of  an  NxM  matrix  P,  and  y,  the  i-th  element  in  a  vector 

T  T  ^ 

Y,  we  can  write  W  =  P  P  and  V  =  P  Y.  The  normal  equations  are  frequently 
written  in  the  form  P^PA  =  p'^Y.) 

The  normal  equations  can  be  solved  for  A  (A  =  W  if  W  is  not 
singular.  (The  BRLESC  routines  for  solving  systems  of  equations  use 
some  form  of  Gauss  elimination. )  The  accuracy  of  this  solution  depends 
on  how  well  conditioned  W  is,  how  large  M  is,  and  the  number  of  digits 
carried  by  the  routine.  The  BRLESC  single  precision  carries  about  l6 
decimal  digits  which  is  equivalent  to  double  precision  on  most 
computers . 

Classical  theory  recognizes  three  situations  for  linear  systems  of 
equations;  (l)  The  equations  are  inconsistent  (e.g.,  x  +  y  =  1, 

X  +  y  =  2),  hence  no  solution  exists.  (2)  The  equations  are  not 
independent  (e.g.,  x  +  y  =  1,  2x  +  2y  =  2),  in  which  case  an  infinity  of 
solutions  exist.  (5)  There  is  a  \mique  solution.  In  theory  it  is 
easy  to  tell  these  cases  apart,  but  in  practice  machine  round-off  error 
blurs  the  differentiation  between  the  three  cases. 

A  frequent  error  in  using  linear  least  squares  is  to  choose  a  set 

of  functions,  q3j^(X),  k=l,2,...,M,  which  are  not  really  independent  over 

the  data  being  fitted.  In  one's  zeal  to  find  an  answer  quickly,  it  is 

easy  to  overlook  even  simple  dependent  relations  among  functions.  The 

particular  set  of  data  and  the  limitation  of  machine  accuracy  sometimes 

make  the  choice  of  clearly  independent  fxmctions  more  difficult.  For 

5  4  4  4 

example,  if  x  is  greater  than  10  ,  x  +  1  =  x  on  BRLESC;  hence,  x  and 

4  5 

X  +1  are  not  independent  on  BRLESC  for  x  greater  than  10  . 

(Similarly,  x  ^  +  1  =  1  for  x  >  10^.) 
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The  easiest  way  to  get  a  linear  least  squares  fit  is  to  use  the 
FORA.ST  LEAST  SQUABES  ERCXjRAM.  This  allows  a  fair  amount  of  flexibility 
with  a  minimum  of  programming  and  gives  "best"  coefficients  and  some 
results  to  check  the  goodness  of  the  fit.  Some  formula  cards  and 
control  cards  must  be  prepared^  but  no  coding  is  needed.  The  MULTIPLE 
REGRESSION  program  is  almost  as  easy  to  use.  It  will  select  a 
statistically  significant  submodel,  from  a  candidate  model  suggested 
by  the  user,  before  carrying  out  the  final  fit.  The  two  subroutines 
GENLSQ  and  POLYLS  (G.L.SQ  and  P.L.SQ  in  FORAST)  are  also  easy  to  use. 

The  user  must  set  up  the  data  in  the  BRLESC  in  a  specified  form,  call 
the  subroutine,  and  then  print,  or  use,  the  results  as  he  wishes.  All 
these  routines  produce  some  resiilts  for  judging  the  goodness  of  the  fit. 
The  FNEQS  and  MATINV  (F.N.E.  and  SY.SNE  in  FORAST)  subroutines  require 
more  programming  but  they  permit  some  additional  freedom. 

2.  Nonlinear  Least  Squares.  If  F(X,A),  the  approximating  function, 

2 

is  not  linear  in  the  a^,  the  minimization  of  U(A)  ==  2ji-l^^i  " 

is  more  difficult.  The  solution  generally  requires  iteration  and  some¬ 
times  the  iteration  fails  to  converge.  Three  of  the  programs  on  our 
list  are  called  "nonlinear  least  squares"  programs.  They  are  quite 
different.  The  simplest  of  the  three  programs  and  the  easiest  to  use 
is  the  FCRAST  NONLINEAR  LEAST  SQUARES  PROGRAM  which  is  similar  to  the 
FQEIAST  linear  LEAST  SQUARES  PROGRAM.  This  program  uses  an  iteration 
method  for  minimizing  U(A)  which  is  sometimes  called  the  Gauss-Newbon 
method,  but  is  locally  called  differential  corrections. 

The  FORTRAN  routine  called  CONSTRAINED  NONLINEAR  LEAST  SQUARES  is 

more  difficult  to  use,  but  it  has  some  features  that  might  make  it 

desirable.  It  permits  constraints  by  adding  the  squares  of  the 

constraints  to  the  sum  of  the  squares  of  the  residuals.  The  resulting 

function  is  minimized  with  a  combination  of  the  method  of  steepest 

19 

descent  and  the  method  of  differential  corrections  that  should  converge 
in  most  cases. 
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The  third  nonlinear  least  squares  routine  is  the  FORTRAN  routine 
called  NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA.  This  routine  uses 
N  data  points  in  the  form:  measured  value 

of  the  i-th  data  point  (Notice  that  no  dependent  variable  is 
prescribed.);  the  standard  error  of  weight  one;  and  R^,  a  KxK 
correlation  matrix  of  cofactors.  The  nonlinear  relation  F(X,A)  selected 
by  the  user  is  assumed  to  be  exact  at  each  of  the  N  points.  The  only 
errors  being  errors  in  each  of  the  NK  variables  (i=l,2, . . . ,N; 

k-l,2,...,K). 

E.  Function  Minimizing  Routines 

If  the  approximating  function  F(X,A)  is  "best"  when  some  function 
U(A)  is  a  minimum^  the  function  minimizing  routines  are  a  rather  natural 
choice  of  routine  for  finding  the  "best"  value  of  A.  In  the  case  of 

least  squares,  U(A)  =  where  R^  is  the  i-th  residual.  This 

choice  of  U(A)  is  frequently  chosen  because  it  is  meaningful  and  easy 
to  analyze.  More  efficient  routines  are  available  for  linear  least 
squares  problems,  but  the  subroutines  FNMIN  and  FDMIN  have  been  used  for 
some  nonlinear  least  squares  problems.  Two  other  fairly  common  choices 

of  U(a)  for  discrete  data  are  U(a)  “^ji-l^^il  U(a)  =  Maximum  |R^|. 

Other  useful  choices  are  possible.  FNMIN  is  the  easiest  of  the  function 

minimizing  routines  to  use.  The  user  must  code  a  subroutine  that 

supplies  U(A)  for  any  given  A.  FDMIN  requires  SU(A)/Sa,,  j=l,2,...,M, 

J 

in  addition  to  U(A).  Neither  of  these  standard  subroutines  permit 
constraints.  One  might  reach  a  satisfactory  solution  for  a  constrained 
problem  by  adding  a  penalty  function  to  U(A),  but  the  nonstandard 
FORTRAN  subroutine  NLEROG  is  probably  a  better  choice.  As  an 
illustration  of  a  penalty  function,  suppose  we  want  to  minimize 

U*(A)  =  ^R?  but  have  constraints  G^(A)  =  0,  r=l,2,...,S.  We  can 

2 

define  U(A)  =  2j1=1^1  Zjr=l'^r^r^'^^  where  the  w^  are  positive  weights. 
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By  adjusting  the  weights  w^^  it  is  possible  to  obtain  an  acceptable 
approximation.  The  CONSTRAINED  NONLUffiAR  LEAST  SQUARES  routine  uses 
such  a  method.  The  NLEROG  routine  uses  this  method  to  handle  equality 


constraints  and  a  sum  of  the  form 


for  inequality  constraints. 


F.  Example 

This  simple  example  was  chosen  to  illustrate  the  effect  of  errors 

in  data  upon  interpolation.  The  table  of  data  (l,l),  (2.1^U), 

2 

(5,9)  was  generated  by  perturbing  the  y  values  of  y  =  x  .  The  data 
could  realistically  be  for  this  function  if  the  readings  were  accurate 
to  within  .5  or  if  they  were  accurate  to  within  .01  with  the  y  values 
for  X  =  2.0  and  x  =  2.1  reversed.  What  is  unrealistic  is  that  the  "true" 
function  is  known. 

The  approximating  functions  corresponding  to  linear,  quadratic,  and 
cubic  interpolation,  the  cubic  spline,  and  the  best  least  squares 
quadratic  were  found.  The  DVDINT  does  not  specifically  display  its 
approximating  function  as  we  do  here. 

Linear  Interpolation 

F(x)  =  3.^x-2.4  0  s  X  s  2 

F(x)  =  -Ux+l2.U  2  X  ^  2.1 

F(x)  =  (50x-69)/9  2.1  s  X  s  3.55 


Quadratic  Interpolation 

F(x)  =  (-7i^0x^+259Ux-17UU0)/ll0  0  ^  X  <  2.1 

F(x)  =  (860x^-3886x+U728)/90  2.1  ^  x  ^  3.55 


Cubic  Interpolation 

F(x)  =  (8o60x^-U7766x^+902UUx-U95l4.8)/990  0  ^  x  ^  3-55 
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Least  Squares  Quadratic 

f(x)  =  .99992l4-44x^-.Oi9998l4-9x+.03854l4-71 


Culaic  Spline  Fit 


f(x)  =  -5. 588964818 (x-l) ^4^. 98896418X-5. 988964818  1  s  x  <  2 

f(x)  =  55.889648i8(x-2.i)^+49.57226005(x-2)^-4.854619082x 

+14.14512781  2  S  X  ^  2.1 

F(x)  =  5. 50802889(3-x)^+10.01705895x-21. 05117688  2.1  s:  x  <  5 

2 

The  following  list  shows  F(x)  and  R(x)  =  f(x)-x  at  x  =  0,  1.5^ 
and  2.5  for  each  of  these  approximating  functions. 


Method 

F(0)  =  R(0) 

F(1.5) 

R(1.5) 

F(2.5) 

R(2.5) 

Linear  Int. 

-  2.4 

2.7 

.45 

6.222 

-  .028 

Quad.  Int. 

-15.8 

4.582 

2.152 

4.511 

-1.959 

Cubic  Int. 

-50.0 

5.605 

5.555 

5.497 

-2.755 

L.  Sq.  Quad. 

.058 

2.258 

.008 

6.258 

-  .012 

Cubic  Spline 

-  2.4 

4.046 

1.796 

4.680 

-1.570 

This  table  shows  that  increasing  the  degree  of  the  interpolating 
polynomial  will  not  necessarily  produce  better  estimates  of  the  "true" 
function.  The  extrapolation  to  x  =  0  demonstrates  why  extrapolation 
from  tabular  data  is  so  •universally  condemned. 

2 

From  the  vie'wpoint  of  approximating  x  from  the  table  of  data,  the 
least  squares  quadratic  was  fairly  successful.  However,  if  we  consider 
this  fit  without  a  priori  knowledge  of  the  "true"  function,  it  too  is 
quite  poor.  The  estimated  standard  error,  ERMS,  is  .572.  This  is  a 
rather  large  error  for  numbers  between  1  and  9*  The  following  table  of 
coefficients  "sigmas  and  t's"  shows  that  the  fit  is  really  not  very 
significant,  and  a^  have  no  significance  and  even  a^  is  of 
doubtful  validity.  (A  frequently  used  rule-of-thumb  is  to  discard  all 
coefficients  with  corresponding  t's  less  than  2.)  The  residuals  for  the 
fo-ur  points  are  -.02,  .40,  -.40,  and  .02,  respectively. 
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j 

a . 

J 

0 . 

J 

t . 

J 

1 

•9999 

.575 

1.738 

2 

-.0200 

2.164 

-.009 

5 

.0383 

2.322 

.018 

With  only  four  points  we  can  tell  very  little  about  the  data  from 
the  residuals,  but  suppose  we  investigated  the  data  and  discovered  that 
the  y^  and  y^  had  been  interchanged.  The  best  least  squares  quadratic 
fit  for  the  four  points  (l,l),  (2,U),  (2. 1,4. 1-),  (5^9)  is 
F(x)  =  1.00496i84x^-.02009923x+. 01536407. 

This  produces  residuals  (.00023,  -.00499,  .00504,  -.00028), 

ERMS  =  .0071,  and 


J 

0 . 

J 

t . 

J 

1 

1.004962 

. 00714 

140.8 

2 

-.020099 

.02881 

-.698 

5 

.01536407 

.02685 

.572 

Again  and  a^  are  not  significant,  but  a^  is.  If  we  refit  without 

a^  and  a^  we  get  F(x)=. 9996245x^  with  residuals  (-.00038,  -.0015,  .0083, 

-.0034)  and  ERMS  =  .0053-  Notice  that  the  residuals  are  generally 

larger,  but  the  unbiased  estimate  of  error  has  dropped  from  .OO7I  to 

.0053  with  the  deletion  of  a„  and  a,. 

2  5 

IV.  DESCRIPTION  OF  ROUTINES 

The  use  of  each  of  the  routines  will  be  discussed  in  this  section. 
The  details  of  Implementing  the  routines  are  presumably  contained  in  the 
program  write  ups,  most  of  which  are  included  in  the  Appendix.  What  I 
hope  to  do  is  briefly  describe  what  the  routine  is  supposed  to  do, 
indicate  the  amount  of  work  necessary  to  use  the  routines,  point  out 
any  unusual  difficulties,  and  add  any  other  remarks  that  seem  pertinent. 
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DVDINT^^^  (D.D.IN^  in  FORAST) 

This  standard  subroutine  uses  Newton's  divided  difference  method  of 
polynomial  interpolation.  The  data  (x^,y^),  i=l,2,...,N,  must  be  ordered 
on  the  dependent  variable  x,  either  strictly  monotone  increasing  or 
strictly  monotone  decreasing.  The  value  of  x  for  which  F(x)j  the 
approximating  function,  is  desired  must  be  between  x^-Cx^-x^)  and 
x^+(x^-x^  .  The  approximating  function  is  an  (M-l)th  degree 

polynomial  through  M  points  near  x.  The  resulting  fiinction  is 
continuous  and  passes  through  each  of  the  data  points.  The  derivative 
however  is  not  continuous.  It  fails  to  exist  at  most  of  the  interior 
points . 


The  most  frequent  error  made  with  this  subroutine  is  attempting  to 
extrapolate  outside  the  acceptable  range.  The  program  produces  an  error 
and  halts  when  this  occurs.  Occasionally  someone  specifies  M  greater 
than  N  and  gets  the  same  error  print.  It  is  my  opinion  that  if  M 
greater  than  three  is  used,  the  user  should  not  only  substantiate  his 
choice,  but  also  produce  a  graph  of  the  fitted  function  with  the  data 
superposed.  The  inspection  of  such  a  graph  can  also  reveal  errors  in 
the  input  data. 


Instructions  for  using  this  subroutine  are  in  the  appendix. 


MATINV,  FNEQS^^^  (S.N.E.,  F.N.E.  etc.^  in  FORAST) 


The  standard  subroutine  for  solving  a  system  of  linear  equations, 
MA-TINV,  can  be  used  to  solve  the  "normal  equations"  associated  with 
linear  least  squares.  For  linear  least  squares  the  approximating 
function  is  linear  in  the  M  unknown  coefficients.  That  is 


F(X,A) 


iaj<P,(X) 


Where  the 


J 


are  appropriately  chosen,  usually 


well  behaved  functions.  We  wish  to  find  A  =  (a^,ag, . . . ,a^)  that 


minimizes  U(A) 


)  .  (y.  -  )  .  ,a  .cp.(X. )) 

Lu=l'''i  AjJ=1  yj'  i" 


Such  an  A  satisfies 
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•  •  « 


L=1  k=l,2,...,M.  These  are  the 


5U(A)/5a^  =  0,  k=l,2,...,M,  or 

,M  „N 

.j=l  -j  L.i=l"k'--i'^j' 
so-called  "normal  equations".  The  formation  of  the  matrix  of 
coefficients  for  these  normal  equations  is  usually  best  accomplished  by 
using  FNEQS,  the  form  normal  equations  subroutine.  These  programs  do 
not  contain  checking  features,  except  for  an  error  print  from  MA.TIWV  if 
a  singular  matrix  occurs.  The  user  is  advised  to  use  GENLSQ  or  POLYLS 
subroutines  or  the  LEAST  SQUARES  PROGRAM  instead. 


One  frequently  made  error  in  using  linear  least  squares  fitting 
procedures  is  to  choose  functional  terms  which  are  not  really  independent 
over  the  given  set  of  data.  The  result  can  be  an  equation  that  fits  the 
particular  set  of  data  used,  but  is  useless  for  intermediate  points.  If 
possible,  it  is  useful  to  reserve  some  data  for  an  independent  check  of 
the  fit.  A  reasonable  rule-of-thumb  is  to  have  three  or  four  data 
points  for  each  \inknown  coefficient.  If  M  is  large,  the  accumulated 
ro-undoff  error  in  inverting  even  well  behaved  matrices  may  produce 
inaccurate  restilts. 


Instructions  for  using  these  subroutines  are  in  the  Appendix. 

GENLSQ  and  FOIYLS  (G.L.SQ  and  F.L.SQ  in  PCRAST) 

The  GENLSQ,  general  least  squares,  routine  is  a  standard  subroutine 

for  making  a  linear  least  squares  fit.  The  input  requirements  are  the 

matrix  (cp.(X.)),  (j=l,2, . . . ,M;  i=l,2,...,N)  and  the  vector  of  the 
J  ^ 

T 

dependent  variable,  (yQ^,y2^  •  •  •  j  where  the  desired  approximating 

function  is  F(X,A)  =  2^^^ajCPj (X)  .  The  POLYLS,  polynomial  least  squares, 
is  a  modification  of  GENLSQ  for  polynomials  in  one  dependent  variable. 

rp 

In  this  case,  the  user  need  only  supply  the  two  vectors  (x^,X2, . . . ,x^) 
and  (y^^^yg, . . .  ,yjj)'^,  and  the  subroutine  automatically  uses 

cpj(Xi)  =  x^  (The  FORAST  version  allows  more  flexibility  in  choosing 
terms  and  permits  weighting  and  omitting  data  points.) 
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These  routines  will  compute  a.,  J=l,2, which  are  "best"  in 

0 


the  least  squares  sense.  That  is,  U(A)  =  2^^_^(y^-F(X^,A) )  i 


IS  a 


minimum.  These  routines  will  compute  the  approximate  values  of  the 
function  F(X.,A),  the  residuals  R.  =  y.-F(X.,A),  the  unbiased  root-mean- 


square  error 


ERMS  =  (  ^^^(R^)/(N-M)f, 


the  estimates  of  error  for  each 


coefficient  a.,  and  an  estimate  of  the  significance  of  the  coefficients 
0 

t .  =  a  ./a . . 

J  0  0 


These  routines  are  easy  to  use  althoxjgh  it  is  easy  to  omit  or 
misplace  one  of  the  many  dummy  variables.  However,  as  with  most 
programs,  the  most  troublesome  and  very  common  errors  are  errors  in  the 
data.  The  problems  with  functional  dependence  and  errors  in  matrix 
inversion  if  M  is  large  are  still  present.  An  independent  check  of  the 
resulting  approximation  is  always  desirable. 

Instructions  for  using  these  subroutines  are  in  the  Appendix. 

FMMIN  and  FDMIN^^^  (FN.MIN^  and  FD.MIH^*^  in  FQRAST) 

These  two  standard  subroutines  are  function  minimizing  routines. 

They  can  be  used  to  attempt  to  solve  function  approximation  problems . 

T 

Suppose  we  wish  to  find  A  =  .  ,b^  which  makes  U(A)  a  minimum. 

These  two  subroutines  attempt  to  minimize  U(A)  directly.  FNMIW  uses 
a  direct  search  technique  which  does  not  require  derivatives.  FDMIH 
uses  a  quasi -Newton  method  which  does  require  first  derivatives  of  U(A) 
with  respect  to  each  of  the  a  . 


The  user  must  supply  the  initial  guess  for  A,  some  control 

parameters,  and  a  subroutine  to  evaluate  U(A),  (U(A)  and  SU(A)/Sa., 

J 

j=l,2,...,M,  for  FDMIN).  The  programs  are  easy  to  use.  However,  they 
are  rather  slow  (compared  to  linear  least  squares)  and  convergence  to 
the  "best"  answer  is  uncertain.  If  the  location  of  the  minimum  of  U(A) 
is  critical,  it  is  desirable  to  check  the  results  of  these  minimizing 


50 


programs.  If  the  function  U(A)  is  sufficiently  well  behaved,  some 
analysis  may  confirm  the  results,  but  this  is  not  always  possible.  One 
partial  check  is  to  repeat  the  minimization  from  several  different 
initial  positions.  If  the  same  solution  is  reached  several  times,  the 
existence  of  at  least  a  local  minimum  there  is  fairly  well  established. 

Most  of  the  difficulties  with  these  two  routines  have  been  from 
errors  in  programming  the  evaluation  of  the  function  or  its  derivatives, 
or  from  using  a  U(A)  which  does  not  have  a  minimum  for  finite  A. 

Neither  of  these  routines  accept  constraints.  Constraints  can 
sometimes  be  included  by  using  a  penalty  function,  but  the  more  general, 
nonstandard,  FORTRAN  subroutine  NLEROG  should  be  a  better  choice. 

Instructions  for  using  these  subroutines  are  in  the  Appendix. 

LEAST  SQUARES  EROGRAM^^,  MULTIPLE  REGRESSIQN^^ ,  NONLINEAR  LEAST  SQUARES^ 

These  three  routines  are  grouped  together  because  of  their  many 
similarities.  All  three  are  complete  FCjRAST  programs  (programmed  by 
L.  W.  Campbell)  which  need  the  same  type  of  formula  definition  cards 
and  control  cards.  All  of  them  have  been  modified  since  their 
respective  descriptions  were  written.  In  particular,  the  restrictions 
on  the  number  of  unknowns  has  been  greatly  increased  in  every  case.  It 
is  desirable  to  get  a  program  deck  directly  from  System  Programming, 
Computer  Service  Division,  along  with  information  as  to  just  what 
options  and  changes  it  contains.  It  might  also  be  wise  to  ask  someone 
from  Systems  Programming  to  check  the  completeness  and  correctness  of 
the  first  set  of  cards  prepared  for  one  of  these  programs . 

Despite  some  uncertainty  as  to  just  what  the  programs  will  do  and 
their  exact  implementation,  these  programs  should  be  given  first 
consideration  for  any  least  squares  fitting  project.  The  fact  that  they 
are  complete  programs  is  both  an  advantage  and  a  disadvantage.  Since 


* 

FORTRAN  programs  for  MULTIPLE  REGRESSION  and  NONLINEAR  LEAST  SQUARES 
are  available. 
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they  are  complete  programs,  no  coding  is  required  (formula  cards  and 
control  cards  are  used  instead)  and  hence,  no  knowledge  of  FQEIAST  (or 
FORTRAN)  is  needed.  Also,  the  results  usually  receive  at  least  a 
cursory  check  before  being  used  for  further  computation.  On  the  other 
hand,  the  resiilts  of  these  programs  are  tabulations  (or  at  best  niombers 
on  cards)  and  are  not  directly  available  for  use  in  other  routines,  as 
the  output  of  a  subroutine  would  be. 

In  all  three  programs,  the  parameters  which  control  the  fitting 
are  punched  on  cards.  Each  data  point,  say  (V1,V2, . . . ,Vk)  is  read  in 
under  the  control  of  a  FORMAT  which  may  have  to  be  changed.  The  user 
defines  the  form  of  the  equation  by  "formula  control"  cards.  These 
formula  control  caxds  may  contain  arithmetic  operations  and  some 
subroutines:  SIN,  COS,  SQRT,  LOG,  EXP,  ARCTAN  and  ARCSIN.  For  example, 
V1*EXP(V6)  =  V2«fr3+L0G(V3)+(VU+V5+SQRT(V7) )$.  This  corresponds  to  the 
formula  Vle^^="a^"+aj^V2^+a2in^(V5)+a^(vU+V5+V7^)  •  (The  constant  term, 
a^,  would  be  automatically  included  in  the  multiple  regression  program, 
■unless  specifically  omitted,  but  not  in  the  other  two  programs.)  The 
output  for  all  of  these  programs  includes  the  coefficients,  approximate 
function  values,  residuals,  ERMS,  and  the  "sigmas  and  t’s"  which  indicate 
the  variance  and  significance  of  the  coefficients . 

The  lEAST  SQUARES  PROGRAM  is  similar  to  the  GENLSQ  subroutine. 

There  is  an  option  which  simplifies  the  input  if  a  polynomial  fit  with 
one  dependent  variable  is  wanted.  The  LEAST  SQUARES  PROGRAM  is  still 
available  but  the  Computer  Support  Division  now  recommends  using  the 
MULTIPLE  REGRESSION  program  instead. 

The  MULTIPLE  REGRESSION  program  treats  the  model  prescribed  on  the 
formula  control  cards  as  a  candidate  model.  The  program  selects  a 
statistically  significant  submodel  and  gives  coefficients,  residuals, 
etc.,  for  this  choice.  A  list  of  the  terms  as  they  are  added  or  removed 
from  the  regression  model  and  the  change  in  ERMS  is  also  given. 
Modifications  for  graphing  and  for  an  input  tape  for  a  companion  program, 
"Prediction  Intervals  for  Estimates  from  Linear  Regression  Models", 
are  described  in  a  description  dated  Nov.  1967- 
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This  very  useful  program  is  thoroughly  described  in  Reference  13 . 

The  NONLINEAR  lEAST  SQUARES  program  is  similar  to  the  lEAST  SQUARES 
PROGRAM.  Since  the  formula  is  not  linear  in  the  unknown  coefficients, 
an  iterative  scheme,  locally  called  "differential  corrections",  is  used 

to  attempt  to  minimize  U(A)  =  2ji=l^yi  ”  is  expanded 

in  Taylor's  series  about  the  latest  estimate  of  A,  say  A,  discarding 
quadratic  and  higher  terms.  AA  =  (Aaj^,Aa2,  •  •  •  ,Aa^)  is  found  so  that 

jSF(X 

linear  least  squares  programs.  Then,  A  is  replaced  by  A  +  jAA, 

(O  ^  p  ^  l)  and  the  process  is  repeated  until  each  j^jj 
prescribed  maximum  nirnber  of  iterations  is  exceeded. 

The  input  for  this  routine  is  more  involved  than  that  for  the  last 

two  programs.  F(X,A),  the  derivatives  SF(X,A)/Sa  ,  j=l,2,...,M,  and  the 

J 

form  y^  -  F  =  ni^st  all  be  described.  The  user  must  also 

supply  (€^,€2^ • • • and  initial  values  of  A.  The  output  of  the  program 
includes  fairly  detailed  results  for  each  iteration  as  well  as  F(X^,A), 
the  residuals  y^  -  F(X^,A),  and  the  "sigma  and  t"  values  for  the  final 
iteration. 

This  program  is  also  fairly  easy  to  use.  The  fact  that  it  is  not 
a  subroutine  is  probably  a  good  thing.  The  resulting  approximation 
should  be  examined  critically  before  it  is  accepted.  If  one  must  have 
a  subroutine,  he  can  code  his  own  program  using  the  differential 
correction  method  and  available  linear  least  squares  subroutines  (or  he 
can  try  one  of  the  function  minimizing  subroutines). 

Instructions  for  using  NONLINEAR  LEAST  SQUARES  are  in  the  Appendix. 
These  instructions  refer  to  Reference  15. 


i,A)/Sa_^)^  is  a  minimuin,  as  in 
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The  author  has  used  the  MULTIPLE  REGRESSION  program,  but  has  no 
direct  experience  with  the  other  two  programs  in  this  section.  Several 
mistakes  were  made  through  carelessness:  A  FORMAT  was  wrong,  a  coznma 
was  omitted  after  a  redefinition  formula,  and  a  plus  sign  was  lost  when 
a  formula  card  was  duplicated.  Despite  these  minor  troubles,  which 
point  out  the  necessity  for  checking  all  computer  input  no  matter  how 
simple,  the  MULTIPLE  REGRESSION  program  was  relatively  easy  to  use. 

PIECEWISE  QUART IC  FIT,  CUBIC  SPLINE^'^,  LEAST  SQUARES  CUBIC  SPLINE 

These  three  routines  have  many  common  features:  They  are  all 
complete  FORAST  programs.  They  are  bivariate  routines  which  serve  as 
alternatives  to  interpolation.  The  approximating  function,  a  piecewise 
low  order  polynomial,  has  continuous  first  derivatives  (the  spline 
functions  also  have  continuous  second  derivatives).  Two  conditions  must 
be  stipulated  for  each  routine  in  addition  to  the  table  of  data  points 
(Xi,yi),  i=l,2,...,N.  The  data  must  be  arranged  so  that  the  independent 
variable,  x,  is  strictly  increasing  (this  may  not  be  strictly  necessary 
for  the  LEAST  SQUARES  CUBIC  SPLINE). 

The  PIECEWISE  QUARTIC  FIT  uses  a  set  of  quartic  equations  to 
approximate  the  table  of  data.  The  derivative  dy/dx  must  be  prescribed 
at  the  first  and  last  point.  The  approximating  function  passes  through 
the  first  and  last  point  and  through  each  "break  point",  the  point  where 
the  approximating  function  changes  from  one  polynomial  to  another.  The 
derivative  at  each  break  point  is  set  equal  to  the  slope  of  the  least 
squares  quadratic  polynomial  involving  the  break  point  and  two  points  on 
each  side  of  it.  The  break  points  are  chosen  so  that  as  many  data 
points  as  possible  are  included  in  the  least  squares  fit  for  each 
succeeding  polynomial  without  exceeding  a  prescribed  root -mean-square 
error.  Finally,  if  the  coefficients  of  the  quartic  equation  are  not 
significant  enough,  the  quartic  polynomial  is  replaced  by  a  cubic 
polynomial.  This  program  should  produce  a  very  useful  fit.  I  have  had 
no  experience  with  it. 


The  CUBIC  SPLINE  produces  an  approximating  function  which  passes 
through  each  data  point,  is  continuous,  and  has  continuous  first  and 
second  derivatives.  Because  of  this  continuity,  this  approximation 
function  may  he  preferable  to  polynomial  interpolation.  However  it  still 
has  the  defect,  in  common  with  polynomial  interpolation,  that  errors  in 
the  data  may  be  magnified  in  the  approximating  function.  The 
approximating  function  is  a  different  cubic  equation  between  each  pair 
of  data  points.  One  obvious  disadvantage  of  this  routine  is  the  creation 
of  a  lot  of  data.  The  table  of  data  containing  2N  numbers  is  replaced 
by  a  table  of  5(N-1)  niunbers:  the  independent  variable  in  the  original 
data,  and  coefficients  for  N-1  cubic  equations. 


The  LEIAST  SQUAEIES  CUBIC  SPLINE  routine  alleviates  the  two  main 
disadvantages  of  the  CUBIC  SPLINE  routine  by  replacing  the  constraint 
that  the  approximating  function  must  pass  through  each  data  point  by  a 
least  square  condition.  The  user  must  decide  how  many  cubic  equations 
he  wishes,  say  K  of  them,  and  define  the  abscissas  of  the  break  points, 
say  x*^,x*2^  •  •  •  It  may  be  difficult  to  choose  good  break  points. 

Common  sense  dictates  that  a  fair  number  of  data  points  be  between  each 
break  point,  that  x*^  ^  x^  and  ^  x^^,  and  that  more  break  points 

are  needed  where  the  slope  of  the  data  is  changing  rapidly.  A  cubic 
spline  approximation,  say  F(x,A),  is  made  through  the  points 
(x*j^,y*j^),  k=l,2, . , .  ,K+1.  The  points  y*^^  are  chosen  by  the  routine  so 


that 


is 


a  minimim. 


The  choice  of  good  break  points  might  prove  to  be  troublesome,  but 
this  appears  to  me  to  be  the  most  attractive  of  the  three  routines 
described  in  this  section.  It  should  be  used  more  often. 

At  the  present  time,  the  only  source  of  these  routines  is  the 
Firing  Tables  Branch  of  EBL.  All  three  programs  are  currently  being 
updated.  Updated  programs  and  descriptions  of  them  may  be  available  by 
the  time  this  report  is  published. 
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POLYGONAL  CURVE 


16 

This  is  a  very  specialized  type  of  cmve  fitting.  The  bivariate 
table  of  data  (x^,y^),  1=1, 2, ...,N  is  replaced  by  a  function,  F(x,A), 
which  is  a  specified  number,  say  K,  of  connected  straight  lines.  The 
POLYGONAL  CURVE  program  chooses  break  points  in  a  very  elegant  manner 
by  combining  dynamic  programming  and  least  squares  to  produce  a 
polygonal  cuive  which  has  the  "best  choice"  of  break  points  (i.e., 

^N  2 

corners)  so  that  2^^_^(y^-F(x^,A) )  is  a  minimum.  Reference  l6  describes 

the  theory  for  these  programs  and  gives  some  examples  of  results. 

Several  FORAST  programs  have  been  coded. 

If  a  program  deck  or  additional  information  is  desired,  contact 
Mr.  C.  M.  Frank,  Army  Materiel  Systems  Analysis  Agency. 

CONSTRAINED  NONLINEAR  LEAST  SQUARES^^ 

This  is  the  only  routine  in  our  list  which  was  not  developed  at 
BRL.  This  routine  is  a  modification  of  a  FORTRAN  Share  program 
developed  by  D.  W.  Marquardt.  It  was  obtained  from  Dr.  Marquardt  by 
Hue  McCoy,  Firing  Tables  Branch,  EBL.  A  few  simple  changes  to  make  the 
program  compatible  with  the  BRLESC  FORTRAN  were  made  by  the  author  of 
this  report.  At  this  time,  the  only  problem  which  has  been  run  locally, 
with  this  program  is  a  simple  test  problem.  This  program  has  been 
included  in  our  list  because  it  has  two  features  not  available  in  o\ir 
other  least  squares  programs:  It  permits  constraints,  and  it  uses  a 
method  of  solution  which  should  be  efficient  and  converge  to  the  "best" 
answer. 

The  constraints  are  included  in  the  problem  by  simply  adding  the 
squares  of  the  unsatisfied  constraints  to  the  usual  least  sqviares 
measure  of  goodness  of  fit.  If  there  are  K  constraints,  G,  (A),  the 

program  minimizes  U(A)  =  where  f^  =  y^-F(X^,A)  for  i=l,2,  ...,N, 

and  ^N+k  ~  ®  ^  satisfied  inequality  constraint. 
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The  method  of  minimizing  U(A)  is  a  combination  of  the  "differential 
correction"  method,  used  with  the  FORAST  NONLINEAR  LEAST  SQUARES  program, 
and  the  well  known  method  of  steepest  descent.  (See  Reference  19*) 


The  user  of  this  FORTRAN  program  must  code  three  subroutines: 

FCODE  to  evaluate  f . ,  PCODE  to  evaluate  Sf./Sa.,  j=l,2,...,M,  and  SUEZ 

1  ^  J 

to  compute  parameters  needed  by  FCODE  and  PCODE. 


This  program  has  several  attractive  features.  The  two  principal 

ones  being  the  ability  to  include  constraints  and  the  higher 

probability  of  convergence  than  the  FORAST  NONLINEAR  lEAST  SQUARES 

program.  The  program  also  supplies  statistical  results  with  which  to 

analyze  the  approximation.  The  user's  subprogram  computes  each  f^  and 

9f. /Sa  .  This  gives  great  flexibility  to  the  definition  of  the 
J 

approximating  function;  hence,  this  program  can  be  very  versatile.  For 
example,  it  should  be  possible  to  use  the  program  to  attempt  to  solve  a 
system  of  nonlinear  equations. 


The  program  has  two  principal  disadvantages:  It  requires  a 
considerable  amovint  of  coding  and  input,  and  anyone  using  the  program  is 
on  their  own  since  no  one  locally  has  had  experience  with  it. 


At  present,  the  only  source  of  a  program  deck  or  description  is 
Hue  McCoy  of  EBL  or  the  author  of  this  report. 

NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA^^ 


This  program  is  \inique  among  the  available  BRIESC  least  squares 

routines  in  that  it  assumes  correlation  of  errors  in  the  data.  (A  basic 

assumption  for  all  the  other  routines  is  that  such  correlation  does  not 

exist.)  This  program  assumes  that  the  approximating  function  F(X,A)  is 

correct  as  to  form,  and  F(X  ,A)  =  0  if  X,  =  (x. ,x,  „, . . .  ,x.„)  and 

1  1  IX  ±d.  IK. 

A  =  (a^,a2, • . . ,a^)  are  correct. 

The  input  for  each  data  point  consists  of  X.  =  (x. ^  ,x.  „, . . .  ,x.,,)  , 
the  i-th  estimate  of  the  variables;  m^,  the  standard  error  of  weight 
one;  and  R^,  the  KxK  correlation  matrix  of  cofactors.  The  routine 
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T 

attempts  to  find  ‘  ^  1=1, 2,..., N  and 

T 

A  =  (a^,a2, . • . ,a^)  so  that  F(X^  +  5^,A)  =  0,  i=l,2,...,N,  and 

U(5)  =  is  a  minimum. 

The  user  must  code  a  MAIN  program  which  supplies  the  data  and 
controls  the  iteration  of  a  subroutine  COLSA  (iteration  is  necessary 
if  F(X,A)  is  nonlinear)  and  the  calling  of  a  subroutine  COLSB  when  a 
satisfactory  fit  is  found.  The  user  must  also  code  a  subroutine  to 
supply  F(X^  +  §^,A)  and  the  first  partial  derivatives  of  this  F  with 
respect  to  k=l,2,...,K,  and  with  respect  to  a.,  j=l,2,...,M. 

It  is  hard  to  draw  any  conclusions  about  a  program  without  using 

it,  but  the  following  things  are  obvious:  A  considerable  amount  of 

coding  must  be  done.  The  program  is  restricted  to  five  independent 

variables  (x. , ,x. . ,x. and  even  with  this  restriction,  51  numbers 
il^  i2^  ^  i5 

are  needed  for  each  data  point.  On  the  other  hand,  the  output  of  COLSA 
and  COLSB  are  extensive  and  fairly  self  explanatory.  Finally,  this  is 
the  only  BRLESC  routine  available  to  perform  least  squares  fits  with 
correlated  data. 

NLPROG^^ 

This  FORTRAN  subroutine  (actually  a  group  of  subroutines)  was 
designed  to  solve  the  general  nonlinear  programming  problem:  Minimize 
U(A)  subject  to  J  inequality  constraints  Gj^(A)  ^  0,  k=l,2,...,J  and  the 
K  equality  constraints  0,  (A)  =  0,  k=J+l, J+2, . . . , J+K.  The  G  (a), 
k=l,2, . . . ,  J+K,  and  U(A)  can  be  nonlinear  functions  of  A  =  (a^,a2,  •  •  •  ,8^^) 
This  problem  is  replaced  by  the  set  of  subproblems:  Minimize  P(A,p)  = 

c  G^(A)/p^  over  A  such  that  Gj^(A)  S  0, 

J  +1 

for  p  =  S  ...  s  until  p^^  is  small  enough. 

There  are  four  versions  of  this  subroutine  depending  on  the  method  used 


58 


to  minimize  P(A,p).  The  user  codes  a  main  program  to  supply  data, 
initial  guesses  at  A,  and  a  number  of  other  parameters.  He  will  also 
have  to  code  a  subroutine  to  evaluate  U{A)  and  the  constraint  equations. 
Depending  on  which  NLTOOG  version  is  used,  the  user  may  have  to  supply 
first  and  second  derivatives  of  U  and  each  of  the  Gj^,  just  first 
derivatives,  or  no  derivatives  at  all. 

This  subroutine  is  nonstandard,  it  may  be  difficult  to  understand, 
it  requires  a  fair  amount  of  programming  and  it  will  be  difficult  to 
check  the  validity  of  the  solution.  On  the  other  hand,  it  can  be 
applied  to  a  wide  variety  of  problems.  This  subroutine  may  converge 
slightly  better  than  the  standard  minimizing  subroutines,  FNMIN  and 
FDMIN,  but  these  two  subroutines  are  much  easier  to  use. 

At  present,  the  only  source  of  program  decks  is  the  author.  A 
complete  description  of  this  routine  is  included  in  Reference  11.  Part 
of  that  description  is  duplicated  in  the  Appendix  of  this  report. 
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APPENDIX:  INSTRUCTIONS  FOR  USING  APPROXIMATION  ROUTINES 


INTRODUCTION 

These  descriptions  are  copied  from  various  published  and  unpublished 
sources.  I  wish  to  thank  L.  W.  Campbell,  CSD,  for  permission  to  use 
material  from  his  publications  and  other  publications  from  Systems 
Programming.  I  also  thank  Dr.  A.  Celnins,  AMD,  for  the  description 
of  his  NONLINEAR  LEAST  SQUARES  ROUTINE  FGR  CORRELATED  DATA,  and 
Joe  Hurff,  EBL,  for  the  description  of  IPAST  SQUARES  CUBIC  SPLINE. 

Standard  FORTRAN  Subroutines 

These  descriptions  are  taken  from  an  impublished  Systems  Programming 
listing  July  25,  I968.  Reference  11  of  this  report  has  a  more  detailed 
description  of  GENLSQ  and  POLYLS. 

LISTING  OF  FORTRAN  SUBPROGRAM  CARD  DECKS  AVAILABLE  FROM 
SYSTEMS  PROGRAMMING, BLDG, 328,  RM  213,  APG,MD. 

KEY  TO  STAT.  NO.  FIELD, 

(l)  INDICATES  INPUT  ARG.,  CALLING  PROG.  SUPPLIES  VALUE 

(R)  INDICATES  RESULT.  SUBPROGRAM  STORES  VALUE  THERE. 

(T)  INDICATES  TEMPORARY  STORAGE. 

(IR)  INDICATES  ARGUMENT  USED  AS  INPUT  AND  RESULT. 

(f)  indicates  ARG.  USED  AS  A  FUNCTION  NAME. 

(S)  INDICATES  ARG.  USED  AS  A  SUBROUTINE  NAME. 

(U)  INDICATES  ARG.  WITH  UNUSUAL  USAGE. 

IMPLIED  TYPE  OF  DUMMY  ARGUMENT  INDICATES  REQUIRED  TYPE  OF  ACTUAL  ARGUMENT, 
EXCEPT  WHERE  NOTED  OTHERWISE. 

IMPORTANT  BRLESCl  RESTRICT  IONS -NO.  OF  DIMENSIONS  MUST  EE  THE  SAME 

BETWEEN  ACTUAL  AND  DUMMY  ARGUMENTS. 

DUMMY  ARRAY  ARG.  CANNOT  HAVE  ACTUAL  ARG. 
THAT  HAS  SUBSCRIPT.  (ACTUAL  ARG.  MUST 
BE  JUST  ARRAY  NAME  WHEN  DUM.  ARG.  IS 
ARRAY.) 

SUBROUTINE  DVDIHT(X,FX,XT,IT,NP,ND) 

C  DOES  DIVIDED  DIFFERENCE  INTERPOLATION. 

C  (l)  X  IS  ARGUMENT  FOR  WHICH  FUNCTIONAL  VALUE  IS  DESIRED. 

C  (R)  FX  IS  NAME  OF  THE  RESULT. 

C  (l)  XT  IS  ARRAY  OF  X  VALUES. (l  DIMENSIONAL) 

C  (l)  FT  IS  ARRAY  OF  FUNCTIONAL  VALUES,  (l  DIMENSIONAL) 

C  (l)  NP  IS  THE  NUMBER  OF  VALUES  IN  XT  AND  FT  ARRAYS. 

C  (l)  ND  IS  THE  NUMBER  OF  POINTS  TO  USE  FOR  EACH  INTERPOLATION. 
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SUBROUTINE  MATINV(A,N,C,NMAX,K,DET) 

C  MATRIX  INVERSION.  A=A**(-l) 

C(IR)  a  is  the  matrix  and  is  REPLACED  BY  ITS  INVERSE. 

C  (l)  N  IS  THE  DIMENSION  OF  THE  MATRIX. 

c  (r)  c  is  used  only  when  k=i  as  described  below. 

C  (l)  NMAX  IS  THE  MAX.  NO.  OF  ROWS  OF  A  AS  DECLARED. 

C  (l)  K  DESCRIBES  OPTIONS.  K=0  MEANS  AN  N  X  N  MATRIX. K=1  MEANS 

C  AN  N  X  N  MATRIX  AND  A  SINGLE  VECTOR  AT  C .  (THE  SOLUTION 

C  VECTOR  REPLACES  THE  C  VECTOR).  K>=2  MEANS  AN  N  X  (N+K-l) 

C  MATRIX.  (THE  (K-l)  VECTORS  ARE  REPLACED  BY  THE  (K-l) 

C  SOLUTION  VECTORS). 

C  (r)  DET  IS  THE  VALUE  OF  THE  MATRIX  DETERMINANT. 

C  WHEN  A  IS  SINGULAR,  AN  ERROR  PRINT  AND  RETURN  WITH  THE 

C  VALUE  OF  DET  SET  TO  ZERO  IS  EXECUTED. 

SUBROUTINE  FNEQS(A,N,C,NMAX,W) 

C  FORM  NORMAL  EQUATIONS  (FULL  N  X  (N+1)  MATRIX). 

C  (R)  a  is  the  matrix  of  normal  EQUATIONS  BEING  FORMED. 

C  A  MUST  BE  CLEARED  TO  ZEROS  BEFORE  FIRST  CALL  OF  FNEQS. 

C  (l)  N  IS  THE  NO.  OF  TERMS  (EXCLUDING  FUNCTIONAL  VALUE). 

C  (l)  C  IS  A  VECTOR  CONTAINING  THE  TERMS  OF  THE  EQUATION  INCLUDING 

C  THE  FUNCTIONAL  VALUE  AS  THE  LAST  TERM. 

C  (l)  NMAX  IS  THE  MAX.  NO.  OF  ROWS  OF  A  AS  DECLARED. 

C  (l)  W  IS  THE  WEIGHT  TO  BE  APPLIED  TO  THIS  EQUATION. 

SUBROUTINE  GENLSQ(X,NRX,F,M,A,NRA,N,C,R,AF,ERMS,SIG,T,DET, IC ) 

C  USES  FNEQS  AND  MATINV  SUBROUTINES. (MUST  INCLUDE  CARDS.) 

C  (l)  X  IS  A  MATRIX  OF  TERMS  OF  EQUATIONS. 

C  (l)  NRX  IS  NUMBER  OF  ROWS  IN  X. 

C  (I)  F  IS  A  VECTOR  OF  FUNCTION  VALUES  FOR  EQUATIONS. 

C  (l)  M  IS  NUMBER  OF  EQUATIONS. 

c  (t)  a  is  a  matrix  of  at  least  (n)x(n+i),is  replaced  by  inverse. 

C  (l)  NRA  IS  NUMBER  OF  ROWS  IN  A. 

C  (l)  N  IS  NUMBER  OF  TERMS  NOT  INCLUDING  FUNCTION  VALUE,  N.LE.99. 

C  (R)  C  IS  A  VECTOR  FOR  N  COEFFICIENTS. 

C  (R)  R  is  a  vector  for  M  RESIDUALS. 

C  (r)  AF  IS  a  vector  for  m  approximate  functions. 

C  (R)  ERMS  IS  THE  ROOT  MEAN  SQUARE  ERRQR,EQUALS  ZERO  IF  M.LE.N. 

C  (R)  SIG  IS  A  VECTCE  FOR  N  SIGMAS. 

C  SIG  IS  INVERSE  ELEMENT  IF  INV.  ELEMENT  IS  NEGATIVE. 

C  (R)  T  IS  A  VECTOR  FOR  N  T  VALUES. 

C  (R)  DET  IS  THE  VALUE  OF  THE  DETERMINANT. 

C  (I)  IC  IS  THE  CONTROL -- 

C  IC  IS  0  COMPUTE  EVERYTHING. 

C  IC  IS  1  COMPUTE  ONLY  COEFFICIENTS. 

C  IC  IS  2  COMPUTE  EVERYTHING  EXCEPT  RESIDUALS  AND  APPROXIMATIONS. 

C  IC  IS  5  COMPUTE  EVERYTHING  EXCEPT  APPROXIMATIONS. 
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SUBROUTINE  POLYLS(X,F,M,A,NRA,N,C,R,AF,ERMS,SIG,T,DET,IC) 

USES  FNEQS  AND  MATINV  SUBROUTINES.  (MUST  INCLUDE  CARDS.) 

X  IS  A  VECTOR  OF  INDEPENDENT  VARIABLE. 

F  IS  A  VECTOR  OF  FUNCTION  VALUES  FOR  POLYNOMIALS. 

M  IS  NUMBER  OF  POLYNOMIALS 

A  IS  THE  MATRIX  OF  AT  LEAST  (N+1)X(N+2)  ,  IS  REPLACED  BY  INVERSE 

NRA  IS  NUMBER  OF  ROWS  IN  A 

N  IS  DEGREE  OF  POLYNOMIAL,  N.LE.  99 

C  IS  VECTOR  FOR  (n+1) COEFFICIENTS 

R  IS  VECTOR  FOR  M  RESIDUALS 

AF  IS  VECTOR  FOR  M  APPROXIMATE  FUNCTIONS 

ERMS  IS  THE  ROOT  MEAN  SQUARE  ERROR, ZERO  IF  M.LE.N+1 

SIG  IS  A  VECTOR  FOR  N+1  SIGMAS. 

SIG  IS  INVERSE  EIEMENT  IF  INV.  EIEMENT  IS  NEGATIVE. 

T  IS  A  VECTOR  FOR  N+1  T  VALUES. 

DET  IS  THE  VALUE  OF  THE  DETERMINANT. 

IC  IS  THE  CONTROL— 

IC  IS  0  COMPUTE  EVERYTHING. 

IC  IS  1  COMPUTE  ONLY  COEFFICIENTS. 

IC  IS  2  COMPUTE  EVERYTHING  EXCEPT  RESIDUALS  AND  APPROXIMATIONS. 

IC  IS  3  COMPUTE  EVERYTHING  EXCEPT  APPROXIMATIONS. 

SUBROUTINE  FNMIN(N,X,FX,FUN,E,EPS,K) 

FINDS  MINIMUM  OF  A  FUNCTION  OF  MORE  THAN  ONE  VARIABLE. 

N  IS  THE  NUMBER  OF  VARIABLES.  N<11  UNLESS  CHANGE  DIMENSION  STATS. 
X  IS  A  LINEAR  ARRAY  CONTAINING  THE  INITIAL  ESTIMATES  OF  THE  N 
VARIABIES  AND  AT  RETURN  CONTAIN  THE  VALUES  AT  THE  MINIMUM. 

FX  IS  WHERE  THE  FUNCTIONAL  VALUE  AT  THE  MINIMUM  WILL  BE  STORED. 

FUN  IS  THE  NAME  OF  A  FUNCTION  OF  2  ARGUMENTS— FUN (X,N)— THAT 
COMPUTES  THE  VALUE  OF  THE  FUNCTION  AT  X.  (AN  EXTERNAL 
STATEMENT  IN  THE  CALLING  PROGRAM  IS  NECESSARY). 

E  IS  THE  NAME  OF  A  SCALAR  WHICH  IS  USED  TO  DEFINE  THE  INITIAL 
TRIAL  STEP  AND  THE  INITIAL  BOUND  FOR  THE  CHANGE  IN  EACH 
VARIABLE.  E>1.  (DELX(i))  INITIAL=E*EPS(I)  AND 
(DELX(I))MAX.  INITIAL=20*(E*EPS(I)). 

EPS  IS  A  LINEAR  ARRAY  OF  N  EPSILONS  DEFINING  THE  ACCURACY 
DESIRED  IN  EACH  OF  THE  VARIABLES. 

K  IF  K=0  INITIALLY,  AN  ERROR  PRINT  AND  HAI/C  WILL  BE 

EXECUTED  WHENEVER  CONVERGENCE  WITHIN  EPS  HAS  NOT  BEEN 
ACHIEVED  AFTER  20+N  ITERATIONS. IF  K  IS  NOT  ZERO  INITIALLY, 
RETURN  IS  EXECUTED  UPON  CONVERGENCE  WITH  K  SET  TO  1, 

OR  AFTER  20*N  ITERATIONS  WITH  K  SET  TO  2. 
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SUBROUTUffi  fdmin(n,x,dx,f,sub,d,eps,epsi,k) 

FINDS  MINIMUM  OF  A  FUNCTION,  USES  DERIVATIVES. 

MUST  BE  FUNCTION  OF  MORE  THAN  ONE  VARIABLE,  I.E.  N.GT.l. 

N  IS  THE  NUMBER  OF  INDEPENDENT  VARIABLES  IN  THE  FUNCTION  TO 

BE  MINIMIZED. N<11  UNLESS  DIMENSION  STATEMENTS  ARE  MODIFIED. 

X  IS  THE  LINEAR  ARRAY  OF  VARIABLES. INITIALLY  CONTAIN  THE 
ESTIMATES  OF  THE  VALUES  AT  THE  MINIMUM. AT  RETURN 
CONTAIN  THE  FINAL  VALUES. 

DX  IS  A  LINEAR  ARRAY  CONTAINING  THE  VALUES  OF  THE  N  PARTIAL 
DERIVATIVES  OF  THE  FUNCTION  EVALUATED  AT  X  BY  THE  SUB 

program.no  initial  values  required. 

F  CONTAINS  THE  VALUE  OF  THE  FUNCTION  AT  RETURN. 

SUB  IS  THE  NAME  OF  A  SUBROUTINE— SUB(N,X,F,DX) —THAT  COMPUTES 
THE  FUNCTIONAL  VALUE  (F)  AND  DERIVATIVES  (DX). 

D  IS  AN  ESTIMATE  OF  THE  IMPROVEMENT  IN  THE  VALUE  OF  THE  FUNCTION 
WHEN  D=0,  ROUTINE  ASSUMES  THE  MIN.  VALUE  IS  NEAR  0. 

EPS  IS  THE  ACCURACY  DESIRED  IN  THE  FUNCTION  VALUE. 

EPSI  IS  A  CONDITION  ON  THE  INDEPENDENT  VARIABLES. 

ABS(DELTAX(i))/ABS(x(i))<EPSI. IGNORED  IF  EPSI  VALUE=0. 

K  IF  K  IS  INITIALLY  ZERO,  AN  ERROR  PRINT  AND  STOP  WILL  BE 

EXECUTED  WHEN  FUNCTION  IS  NOT  CONVERGING. IF  K  IS  NOT  ZERO 
INITIALLY, RETURN  IS  EXECUTED  WITH  K  SET  TO  1  WHEN 
CONVERGENCE  IS  SATISFIED  OR  K  SET  TO  2  WHEN  THERE  IS  NOT 
CONVERGENCE. 


Standeird  FORAST  Subroutines 


Most  of  these  descriptions  are  taken  from  BRL  Report  No.  1275, 

The  FORAST  Programming  Language  for  ORDVAC  and  BRLBSC  (Revised) ,  by 
L.  W.  Campbell  and  G.  A.  Beck,  March  1965*  The  descriptions  of  FN.MIN 
and  FD.MIN  are  unpublished  documents  available  from  Systems  Programming. 


D.D.IN)X)FX)Xo)Fo)tpt)n)ix) 

X 

FX 

(Divided  Difference  Inter¬ 

Xo 

polation) 

Must  use  all  three 

Fo 

optional  arguments  or 

none.  If  omitted. 

tpt 

(5)1)1)  is  used. 

D.D.SX)Fo)FX)^ 

n 

Use  this  to  interpolate 

more  functions  using  the 

ix 

is  the  address  of  the  argument, 
is  the  address  of  the  result, 
is  the  initial  address  of  the 
table  of  Xi ' s . 

is  the  initial  address  of  the 

table  of  Fi ' s . 

is  the  number  of  entries  in 

the  table,  (no.  of  Xi's) 

is  the  number  of  points  to  use 

in  the  interpolation. 

is  the  distance  between  entries 


same  value  of  X.  (must 
use  enter) . 


in  the  X  table 


if  is  the  distance  between  entries 
in  the  F  table. 

The  following  matrix  manipulation  routines  are  available: 

S.N.E  )Al,l)n)Co)DET‘^  Al,l:  Bl,l:  Cl,l  are  addresses 

MAT.INV.  )Al,l)n)Co)DET^  of  the  first  elements  of  matrices. 

To  omit  the  use  of  n  is  the  number  of  unknowns  (rows). 

Co  and  use  DET  it  is  Co  is  the  address  of  the  first 

necessary  to  write  element  of  the  solution. 

Al,l)n))DET)  DET  is  the  address  of  the  determinant. 

SY.SNE  )Al,l)n)Co)DET'^  Cl  is  the  address  of  the  first 

SY.INV  )Al,l)n)Co)DET^  CAP  coefficient  of  the  given  equation. 

F.W.E.  )Al,l)n)Cl)^  W  is  the  address  at  weights  for- the 

F.O.MAT  )Aljl)n)^  equation. 

MAT.M  )Al,l)Bl,l)Cl,l)i)  i  is  the  number  of  rows  in  A(or  AT). 

j)k)^  j  is  the  number  of  cols,  in 

A (or  AT)  and  is  equal  to  the 
no.  of  rows  in  B(or  bT) . 
k  is  the  nmber  of  columns  in 
B(or  bT). 

Additional  comments  on  the  above  matrix  subroutines:  The  S.N.E. 
(Solve  normal  equations)  assumes  all  elements  of  a  matrix  having  n  rows 
and  n  +  1  columns  are  stored  in  the  memory  by  rows.  The  SX.SNE 
(symmetric  solve  normal  equations)  assumes  that  only  the  upper  triangle 
of  an  n  X  n  +  1  matrix  is  stored  and  SX.INV  (symmetric  inversion) 
assumes  that  only  the  upper  triangle  of  an  n  x  n  matrix  is  stored. 
S.N.E.;  MAT.INV;  SY.SNE;  and  SY.INV  all  replace  the  original  matrix 
with  its  inverse.  The  SY.SNE  stores  the  solution  vector  only  at  Co. 

The  F.N.E.  (form  normal  equations)  assumes  that  the  upper  triangular 
augmented  matrix  has  been  cleared  by  the  program  before  it  is  entered 
with  the  first  equation.  The  F.N.E.  produces  a  matrix  that  can  be 
solved  with  the  SY.SNE.  The  F.O.  MAT  (fill  out  matrix)  will  take  an 
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augmented  upper  triangular  matrix  (as  generated  by  F.N.E.)  and  replace 
it  with  an  augmented  square  matrix  (as  needed  by  S.N.E.)* 

The  S.N.E.  will  attempt  to  rearrange  rows  of  the  matrix  when  it 
finds  a  zero  diagonal  element  while  it  is  computing  the  inverse.  The 
row  rearrangement  does  not  affect  the  arrangement  of  the  solution  vector, 
however  the  inverse  matrix  will  not  be  correct  if  any  rows  were  actually 
rearranged.  Rearrangement  can  be  avoided  by  use  of  the  "not"  option 
as  explained  below. 

Additional  BRLESC  S.N.E.  options: 

S.N.E. )Al.l)n)Co)DET)drow)dcol)Bl)db)dc)ZERO)not)^ 

If  "drow"  is  specified,  it  is  the  spacing  between  rows;  i.e.  the 
address  A2,l  -  address  Al,l. 

If  "dcol"  is  specified,  it  is  the  spacing  between  columns  (which 
is  the  same  as  spacing  between  elements  within  a  row) . 

If  B1  is  specified,  the  n  positions  beginning  at  B1  are  used  as 
the  column  vector  instead  of  the  (n+l)  column  of  the  matrix. 

If  "db"  is  specified,  it  is  the  spacing  between  the  elements  of 
the  coltmm  vector. 

If  "dc"  is  specified,  it  is  the  spacing  between  elements  of  the 
solution  vector. 

If  ZERO  is  specified,  it  is  the  address  of  the  number  which  will  be 
used  to  check  for  zero  diagonal  elements.  Those  diagonal  elements  whose 
absolute  value  are  less  than  ZERO  will  be  considered  as  zero  for  the 
rearrangement  test. 

If  "not"  is  any  address  different  from  zero,  the  S.N.E.  will  not 
rearrange  any  rows. 

When  any  or  all  of  these  spacing  options  are  omitted  (or  zero), 
the  normal  consecutive  spacing  of  elements  is  assumed. 
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For  MT.INV  "dxow",  "dcol",  ZERO,  and  "not"  may  be  specified 
when  needed  and  have  the  same  meaning  as  for  the  S.N.E.  except  "not" 
has  the  opposite  meaning.  MT.INV  does  not  normally  rearrange  any 
rosw  and  will  do  so  only  when  "not"  is  specified  as  non-zero. 

Note  that  when  optional  axldresses  are  omitted  any  place  except  at 
the  end  of  an  ENTER  statement,  the  right  parenthesis  must  still  be 
written  for  each  omitted  address.  In  particular,  the  above  options  for 
the  MAT.INV  subroutine  must  correspond  to  the  same  position  on  the  list 
of  addresses  as  used  by  the  S.N.E.  since  they  are  just  different 
entrance  points  to  the  same  subroutine. 

G.L.SQ 

or 

P.L.9a)X)ix)F)if)m)Al,l)n)C)R)ir).AP)iaf)ERMS)SIG)T)DET)w)iw)EQSEQ)TSEQ  $ 

(General  or  polynomial  least  squares  data  fitting. ) 

X  For  G.L.SQ,  X  is  the  location  of  the  first  term  of  the  first 
equation.  Terms  must  be  stored  consecutively. 

For  P.L.SQ,  X  is  the  first  independent  variable. 

ix  For  G.L.SQ,  ix  is  the  distance  from  one  equation  to  the  next 
one.  For  P.L.SQ,  ix  is  the  distance  from  one  independent 
veiriable  X  to  the  next  one. 

F  is  the  function  value  for  the  first  equation  or  polynomial. 

if  is  the  distance  between  function  values. 

m  is  the  actual  total  number  of  equations  or  "points"  that  are 
to  be  used  in  computing  the  fit.  (it  must  not  Include  those 
skipped  by  using  EQSEQ.) 

Al,l  is  a  block  of  storage  that  must  be  large  enough  for  an 
augmented  (n  x  n)  symmetric  matrix. 

n  For  G.L.SQ,  n  is  the  actual  number  of  terms  to  be  used  in 

each  equation,  (it  must  not  include  those  skipped  by  using 
TSEQ.)  For  P.L.SQ,  n  is  one  less  than  the  number  of  terms 
and  is  the  degree  of  the  polynomial  when  all  the  terms  are 
used. 

C  is  the  initial  address  for  consecutively  storing  the  n 

coefficients,  (if  n  ^  58,  n  +  1  spaces  must  be  allowed  at  C.) 
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R  is  the  initial  address  for  storing  the  m  residuals. 

ir  is  the  distance  desired  between  residuals^  i.e.  the 

increment  for  the  R  address. 

AT  is  the  initial  address  for  storing  the  m  approximate  function 
values . 

iaf  is  the  increment  for  the  AF  address. 

ERMS  is  the  store  address  for  the  root -me an -square  error.  Zero  is 

stored  when  m  ^  n  or  when  ^  ^  n. 

SIG  is  the  initial  address  for  consecutively  storing  the  n  "sigmas". 

SIG.  =  ERMS  *  SQRT(inv.el.A.  .) 

(if  the  inverse  element  ^  is  negative,  it  is  stored  for 

SIG.  and  T.  =  0. ) 

1  1 

T  is  the  initial  address  for  consecutively  storing  the  n  "t's". 

T.  =  C./SIG. 

Ill 

DEI  is  the  address  to  store  the  determinant. 

W  is  the  initial  address  of  the  weights  to  be  used. 

iw  is  the  increment  for  the  W  address. 

EQSBQ  is  the  initial  address  of  a  consecutive  sequence  of  numbers 
that  have  a  one  to  one  correspondence  with  each  equation 
(or  point)  stored  at  X.  A  zero  number  indicates  that  the 
corresponding  equation  (or  point)  is  to  be  used  and  a  non-zero 
number  indicates  that  it  should  not  be  used.  Note  that 
this  sequence,  if  used,  must  contain  m  zero  numbers. 

TSEQ  is  the  initial  ad.dress  of  a  consecutive  sequence  of  numbers 
that  have  a  one  to  one  correspondence  with  the  terms  in  each 
general  equation  or  with  the  powers  of  X  in  a  polynomial.  A 
zero  number  indicates  that  the  corresponding  term  or  power 
of  X  should  be  used  and  a  non-zero  number  indicates  that  it 
should  not  be  used.  Note  that  this  sequence,  if  used, 
must  contain  n  zero  numbers. 
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A  FUNCTION  MINIMIZING  SUBROUTINE  WITHOUT  CALCULATING  DERIVATIVES 


Donald  Taylor  and  John  Wortman  Jan  1966 

The  subroutine  will  find  the  minimum  of  a  function  f(X,  ....X  ).  It 

1  n 

requires  that  the  programmer  provide  for  the  evaluations  of  the  function. 
(At  present,  available  on  BRLESC  only). 

The  entrance  sequence  is  as  follows: 

ENTER(FN.MIN)n)X^)F)EVF)E)H^)EPS^)ERR) 


n  is  the  number  of  independent  variables  in  the  function  to 
be  minimized. 

X^  is  the  address  of  the  first  of  n  consecutive  positions  for 
the  independent  variables.  Initially,  these  positions 
should  contain  estimates  of  the  values  at  the  minimima; 
upon  exit  from  the  subroutine,  they  contain  the  final 
values  of  the  independent  variables. 

F  is  the  address  of  the  value  of  the  function. 

EVF  is  the  address  of  the  Programmer' s  function  evaluation. 

It  must  use  X^....X^  as  the  values  of  the  independent 

variables  and  store  the  resulting  function  value  at  F. 

Use  the  name  (FN.RET)  to  return  to  the  subroutine. 

E  is  the  adLdress  of  a  number  which  is  used  to  define  the 
initial  trial  step  and  the  initial  bound  for  the  change 
in  each  independent  variable.  Normally,  E  >  1. 

initial  =  E  x  EPS^  (^^)max.  initial  =  20  (E  x  EPS^) 

2 

is  the  address  of  the  first  of  n  +  4n  consecutive 
positions  used  for  temporary  storage.  (Direction 
matrix,  etc.) 

EPS^  is  the  address  of  the  first  of  n  consecutive  positions 
which  should  contain  the  accuracy  desired  in  each  of  the 
independent  variables. 

ERR  OPTIONAL.  The  subroutine  will  send  control  to  this 

address  (instead  of  producing  an  ERROR  print)  if,  after 
20n  iterations,  the  requested  acc-uracy  has  not  been 
satisfied. 
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It  is  possible  to  maximize  a  function  f  by  minimizing  the  function 
-f.  In  least  squares  fitting,  the  subroutine  has  been  used 
successfully  to  fit  some  non-linear  functions  for  which  'differential 
corrections'  did  not  work  well.  In  this  case,  the  function  to  minimize 
is  the  sum  of  the  squares  of  the  residuals. 


The  subroutine  is  derived  from  a  method  described  by  M.J.D. Powell. 
It  is  a  variation  of  the  well  known  method  of  minimizing  a  function  of 
several  variables  by  changing  one  parameter  at  a  time. 

The  method  does  not  recognize  constraints  on  the  variables. 
Sometimes,  it  may  be  possible  to  apply  some  constraints  within  the 
function  evaluation  program  by  assigning  some  relatively  large  value  to 
the  function  whenever  the  constraint  has  been  violated.  There  is, 
however,  no  assurance  that  this  will  succeed. 

The  error  print  when  the  number  of  iterations  exceeds  20  n  (unless 
bypassed  by  the  option)  is  ITER.>20xN  and  the  number  printed  is  the 
minimum  value  of  the  function  at  this  point. 

It  is  sijggested  that  several  sets  of  initial  values  of  the 
independent  variables  should  be  tried  to  see  that  they  converge  to  the 
same  minimum  to  give  some  assurance  to  the  result.  The  subroutine  does 
not  do  this. 


A  FUNCTION  MINIMIZING  SUBROUTINE  (WITH  DERIVATIVES) 

J.C.  Torrey  and  John  Wortman  November  I965 


This  subroutine  will  find  the  minimum  of  a  function  f(x-,...x  ). 

^  1  n' 

The  routine  is  called  ED. MIN  to  emphasize  that  the  programmer  must 
provide  evaluations  for  the  partial  derivatives  of  his  function  as  well 
as  for  the  function  value.  At  present  it  is  available  only  on  ERLESC. 

The  entrance  sequence  is  as  follows: 

ENTER(FD.MIN)n)Xl)EVF)D)EPS)Hl)F)EFSl)^  ,  where 

n  is  the  number  of  independent  variables  in  the  function 
to  be  minimized. 


M.J.D. Powell;  An  Efficient  Method  for  Finding  the  Minimum  of  a  Function 
of  Several  Variables  Without  Calculating  Derivatives.  The  Computer 
Journal,  Vol.  No.  2,  July,  1964. 
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XI  is  the  address  of  the  first  of  n  consecutive  positions 

for  the  independent  variables.  Initially,  these  positions 
should  contain  the  programmer's  estimates  of  the  values 
at  the  minimum;  when  the  routine  is  finished,  they  hold 
the  final  values  of  the  independent  variables. 

EVF  is  the  address  of  the  programmer's  function  and  derivative 
evaluation.  Retvirn  from  EVF  to  the  routine  is  done  by 
GOTO (fD. RET )^. 

D  is  the  address  of  an  estimate  of  the  improvement  in  the 
value  of  the  function.  This  estimate  will  not  be 
critical  except  in  cases  where  the  routine  is  to 
distinguish  between  local  minima.  (Option:  When  D  =  0, 
the  routine  assumes  that  the  function' s  minimum  value  is 
near  zero  and  ^  zero.) 

EPS  is  the  address  of  the  accxiracy  desired  in  the  fvmction 
value . 

2 

HI  is  the  address  of  the  first  of  n  +  9n  consecutive 
positions,  used  for  temporary  storage. 

F  is  the  address  of  the  first  of  n+1  positions  for  the 
value  of  the  function  and  its  n  partial  derivatives. 

(Note  that  the  function  is  stored  at  F,  with  derivatives 
following  in  the  same  order  as  the  corresponding 
variables  at  XI.) 

EPSl  (optional)  is  the  address  of  a  cauchy-like  condition  on 

the  independent  variables.  At  the  last  routine  iteration, 

<  EPSl,  i=l...n 


The  programmer  will  use  the  EPSl  option  when  his  interest  is  in  the 
values  of  his  variables  at  the  minimum,  rather  than  the  function  Itself. 

Note  that  FD.MIN  can  be  used  to  maximize  a  function.  To  maximize 
f,  minimize  -f. 

The  routine  has  been  used  successfully  in  least  squares  fitting 
that  would  normally  be  done  by  differential  corrections.  To  do  this, 
the  function  to  minimize  is  the  sum  of  the  squares  of  the  residuals. 

The  value  of  D  at  entrance  should  be  zero. 


FD.MIN,  derived  from  a  routine  described  by  Davidon  ;  uses  a 
variable  n  x  n  matrix  M  as  a  metric  in  searching  for  the  best  values  of 
the  X^.  The  programmer  desiring  to  impose  constraints  on  his  variables, 

or  to  use  information  about  them  to  speed  the  search,  should  modify  M 
the  first  (but  only  the  first)  time  the  routine  enters  his  function 

2 

evaluation.  M  is  stored  in  the  first  n  positions  of  the  block  HI,  and 
is  set  to  the  linit  matrix  at  the  start  by  the  subroutine. 

Setting  the  diagonal  elements  of  M  to  the  squares  of  the  estimated 
error  in  the  initial  values  of  the  variables  may  effect  a  substantial 
increase  in  the  speed  of  minimization.  Thus,  if  =  l4  ±  4,  then 

m^  ^  =  l6j  but  if  X^  =  l4  ±  .1,  then  m^  ^  =  .01. 

If  the  programmer  desires  to  impose  linear  constraints  on  his 
variables,  he  modifies  the  matrix  at  the  first  function  evaluation. 

For  the  constraints 


1  =  “i ' 

^  b^Xi  =  kg  ,  etc . , 

he  must  choose  M  so  that 

Za.m.  .  =  0 
1 

I  ^i“^ij  "  °  ^ 

and  initial  values  of  the  X^  to  satisfy  the  constraints. 

The  subroutine  has  an  error  exit  which  prints  'FDMIN  DOWN'  for  its 
error  number.  A  faulty  derivative  evaluation  is  the  most  likely  cause 
of  error,  but  the  programmer  should  assure  himself  that  his  fiinction  has 
a  minimum. 


Finding  a  minimm  of  a  general  function  is  an  -uncertain  process. 
Davidon  suggests  converging  to  the  minimum  several  times  from  varying 
initial  values  before  accepting  it.  FD.MIN  does  not  do  this,  but  its 
users  may  add  assirrance  to  their  results  by  following  his  lead. 


William  C.  Davidon,  "Variable  Metric  Method  for  Minimization", 
ANL-5990,  Physics  and  Mathematics  (TID-4500,  l4th  ed.)  AEC  R&D  Report. 
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FQRAST  Programs 

The  LEIAST  SQUARES  PROGRAM  has  been  extended  considerably  since  its 
original  description  and  amendment  (Reference  12)  were  written  in  1965- 
This  program  is  still  available,  but  the  Computer  Support  Division 
recommends  using  the  MULTIPLE  REGRESSION  program  instead. 

The  description  of  the  MULTIPLE  REGRESSION  program  (Reference  15) 
is  much  to  bulky  to  include  here.  The  description  of  the  NONLINEAR 
LEIAST  SQUARES  PROGRAM  is  included.  This  description  refers  to  Reference 

15. 

FORTRAN  programs  and  descriptions  for  both  MULTIPLE  REGRESSION  and 
NONLINEAR  LEAST  SQUARES  will  be  available  by  the  time  this  report  is 
published. 

The  LEIAST  SQUARES  CUBIC  SPLINE  is  the  only  other  description 
available  at  this  time  for  programs  in  this  section. 

The  PIECEWISE  QUARTIC  FIT,  and  the  CUBIC  SPLINE  have  been  used  by 
various  members  of  the  Firing  Tables  Branch  of  EBL.  They  may  now  have 
programs  and  instructions  for  them  as  they  have  for  the  LEAST  SQUARES 
CUBIC  SPLINE. 


NONLINEAR  LEAST  SQUARES  PROGRAM 

L.  Campbell  December  I966 


The  multiple  regression  computer  program  (as  described  in  BRL 
Report  No.  1550)  has  been  used  as  a  basis  for  developing  a  new  general 
purpose  non-linear  least  squares  program  for  BREESC.  Except  for  the 
exceptions  noted  here,  the  rules  for  using  this  program  are  the  same  as 
for  the  regression  program. 

For  the  non-linear  least  squares,  the  expression  that  is  to  be 
used  to  fit  the  observed  data  must  be  included  as  a  redefinition 
formula  except  it  must  be  called  F.  The  right  side  of  this  formula 
must  be  inclosed  in  parentheses  if  it  has  +  signs  between  "terms". 
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For  example: 


F  =  (A  +  EXP(A1*V1))  , 

where  A,A1,...,A24  must  be  used  to  refer  to  the  coefficients  that  are  to 
be  determined  from  initial  approximations  to  these  coefficients.  (A  and 
Ao  both  refer  to  the  first  coefficient.) 

SF  9F 

The  partial  derivatives  ^  , . must  also  be  defined  by 

redefinition  formulas  and  they  may  be  called  P,P1,P2,...  (or  P's  may  be 
used  if  R's  are  also  used  in  the  main  formula).  For  example,  the  partial 
derivatives  of  the  above  F  expression  would  be  written  as: 

P  =  1,  PI  =  VI  *  EXP(A1*V1)  , 

The  main  formula  must  define  a  residual  as  the  sum  of  all  the 
partial  derivatives,  so  on  the  left  of  the  =  symbol  should  be  "some 
V(or  R)  -  F"  and  on  the  right,  each  i-th  term  is  the  name  used  for  the 
partial  derivative  with  respect  to  the  i-th  coefficient.  The  main 
formula  for  the  above  example  and  using  V4  as  the  observed  value  of  the 
function  would  be 


Vit-  -  F  =  P  +  Pl^ 

Following  the  main  formula  must  be  one  or  more  lines  of  initial 
estimates  for  the  coefficients  punched  in  10  column  fields  with  eight 
per  line.  They  should  be  punched  with  decimal  points  and  may  have 
exponents. 

Following  the  coefficient  estimates  must  be  epsilons  for  testing 
convergence  of  the  coefficients.  Each  coefficient  must  have  its  own 
epsilon  and  they  are  to  be  punched  in  10  column  fields  with  eight 
per  line  and  may  have  decimal  points  and  exponents. 

The  data,  with  optional  header  lines,  follows  and  it  must  be 
followed  with  2  blank  lines.  Data  may  be  on  tape  or  cards  or  previous 
data  may  be  used  by  using  just  2  blank  lines,  just  like  the  regression 
program.  The  initial  value  of  AH  is  1  in  this  program  too,  so  a  single 
header  card  is  expected  to  precede  the  data.  The  amovint  of  data 
allowed  on  tape  is  unlimited,  the  amount  of  data  allowed  on  cards  is 
limited  to  the  BRLESC  memory  capacity  -  9000  approximately. 

The  program  limits  the  number  of  iterations  to  20  unless  an 
"mi  =  i"  control  card  is  used  where  i  is  a  new  maximum  number  of 
iterations.  When  convergence  is  not  reached  within  this  maximum 
number  of  Iterations,  the  program  prints  "DIDN'T  CONVERGE  IN  i 
ITERATIONS."  and  will  not  print  any  residuals,  sigmas  or  t's. 
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An  "USE  DELTA*  p"  control  card  may  be  used  to  cause  each  change  in 
the  coefficients  to  be  less  (or  more)  than  what  the  program  computes. 
Using  p  <  1  may  allow  convergence  when  p  =  1  (its  normal  value)  will 
cause  divergence,  (p  >  1  should  not  be  used.)  Each  computed  change  in 
each  coefficient  is  multiplied  by  p  before  it  is  used  to  actually  change 
the  coefficient  for  the  next  iteration. 

There  is  no  rescaling  in  this  non-linear  program.  All  residuals 
and  approximate  functions  printed  will  be  in  their  redefined  scale. 

The  output  Includes  for  each  iteration,  the  current  root-mean- 
square  error,  the  maximum  residual,  the  coefficients  before  being 
changed,  the  "delta"  or  amount  of  change  for  each  coefficient  and  the 
new  coefficients. 

After  convergence,  the  program  normally  prints  all  the  approximate 
values  of  the  function  and  residuals,  the  same  as  the  regression  program 
except  there  is  no  rescaling.  The  final  root -mean-square  error  and 
maximum  residual  is  printed  and  the  "sigmas"  and  "t's". 

The  following  control  cards  may  be  used  and  have  the  same  meaning 
as  in  the  regression  program. 

PR.RES.>= 

NO  RESIDUALS 
RESIDUALS 
STOP  ERMS= 

ANGLES  ARE  IN  MILS 
ANGLES  ARE  IN  DEGREES 
ANGIES  ARE  IN  RADIANS 
VW= 

VN= 

V.= 

IDE 
AH= 

New  Control  Cards: 

MI=  i 


USE  DELTA*  p 


where  i  is  maximm  no.  of  iterations 
allowed,  (i  =  20  initially.) 

where  p  is  multiplied  times  actual  computed 
coefficient  "deltas"  before  actually  computing 
new  coefficients. 


CARD  INPUT 
TAPE  INPUT 
TAPE  MF  TO 
TAPE  MFMF 
SAME  FORMULAS 
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Input  Sequence; 


Control  Cards 

Re-definition  formulas  (must  include  F  =  formula) 
Main  formula  (Observed  Function  -  F  =  Partial 
derivatives . ) 


or  SAME  FORMULAS 


Initial  Coefficient  Estimates  (lO  col\mm  fields) 


Epsilons  for  convergence  for  each  coeff.  (lO  column  fields.) 

Heater  cards,  if  AR  ^  0  Omit  to  re-use  same  data. 

Data  Data  may  be  on  tape  or  cards 

2  blanks 

(May  repeat  above  sequence  any  number  of  times.) 


Output  Sequence: 


Header  lines,  if  any. 

First  12  data  numbers . 

Two  lines  of  program  parameters. 
All  formulas, 
blank  line 


Current  ERMS  and  Max.  Residual 
Initial  Coefficients 

Computed  changes  of  coefficients  repeated  for  each 

If  p  7^  1,  actual  changes  of  coefficients  iteration. 

New  coefficients 
blank  line 


Approximate  function  and  Residuals 
or 

Data  line  no.,  original  function,  and  residual. 
Final  ERMS  and  final  Max.  Residual 


Sigmas 

t's  (If  t  =  0,  then  corresponding  sigma  is  negative  diag.  matrix 
element. ) 

(Next  output  sequence  will  start  on  a  new  page.) 
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PROGRAM  CUBIC  SPLINE  LEAST  SQUARE 

PROGRAMMER  JOE  HURFF 

DATE  MAY  I969 

LANGUAGE  FORAST 

DESCRIPTION  FOR  GIVEN  DATA  SETS,  THIS  IROGRAM  YIELDS  AN  APPROXIMATING 

FUNCTION  IN  THE  FORM  OF  A  SERIES  OF  CUBIC  POLYNOMIALS. 

THE  FUNCTIONS  ARE  DERIVED  BY  USING  BOTH  THE  PRINCIPLE  OF 
LINEAR  LEAST  SQUARES  AND  CUBIC  SPLINE  FITTING.  THIS 
PROCESS  MAY  USUALLY  BE  USED  SUCCESSFULLY  FOR  SMOOTHING 
AND  PROVIDING  FIRST  AND  SECOND  DERIVATIVES. 


PROCEDURE 


I  DATA  CARDS  (l2  DIGIT  FLOATING  POINT)  THERE  MAY  BE  A 
MAXIMUM  OF  500  DATA  POINTS 

1.  CARD  SET-UP 
COLS 

1-12  Xj 

13-2L  F  (Xj) 

25-80  BLANK 


If  there  are  fewer  than  10  break  points,  the  S  must  follow  in  field 
immediately  after  the  last  one. 
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II  HEADER  CARDS  (NOIffi) 


TIME 

OUTPUT 


III  CONTROL  CARD  (l2  DIGIT  FLOATING  POINT) 

COLS 

1-12  F'(Xq) 

13-24  F'(Xjj) 

73-75  3  DIGIT  CODE  -  USER'S  OPTION 

{0  -  FIT  F’  (Xq)  exactly 

1  -  CONSTANT  F”  (X)  OVER  FIRST  INTERVAL 
2  -  FIT  F(Xq)  EXACTIY 


77 


b  -  FIT  F' (X^)  EXACTLY 

^1  -  CONSTANT  F"  (X)  OVER  THE  LAST  INTERVAL 
2  -  FIT  F(X^)  EXACTLY 


78-80  3  DIGIT  CODE  -  USER'S  OPTION 


IV  ORDER  OF  INKTT 

1.  CONTROL  CARD 

2.  BREAK  POINT  CARDS 

3.  DATA  CARDS 

4.  BLANK 

5.  PROB  CARD 

APPROXIMATEIY  ONE  HALF  MINUTE  PER  CASE 


THE  OUTPUT  IS  COMPLETEIY  LABELED 


Other  FORTRAN  Routines 

No  description  of  the  Share  program  ve  called  CONSTRAINED  NONLINEAR 
LEAST  SQUARES  Is  given.  Dr.  Celmlns  kindly  prepared  a  description  of 
the  NONLINEAR  LEAST  SQUARES  FOR  CORRELATED  DATA  routine.  The 
description  of  NLPROG  was  taken  from  BRL  Memorandum  report  No.  1958. 
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NON-LINEAR  LEAST  SQUARES  SUBROUTINES  FOR 
GCRRELATED  DATA  (COLSA,  COLSB) 

A.  Celmins  December  I968 

Reference:  (I8)  BRL  Report  I36O  (March  I967),  Appendix 


Problem  Outline 

A  general  least  squares  problem  can  be  formulated  as  follows:  We 
have  observed  r  sets  of  n  quantities  Between  these 

quantities  a  functional  relationship  which  depends  on  m  imknown 
parameters  . . . ,y^  is  assumed: 


n  i  m 


(1) 


If  r  >  m,  we  compute  the  parameters  and  some  corrections  of  the 
observations  such  that  Equation  (l)  is  satisfied  at  all  r  observation 
sets  and  the  sum  of  correction  squares  is  a  minimum.  For  this 
computation  the  correction  squares  should  be  weighted  differently 
depending  on  their  accuracies  and  on  correlations  between  them. 


The  subroutine  COLSA  computes  the  values  of  the  parameters  y., 

J 

their  accuracies,  and  their  cofactor  matrix  (i.e.  variance -co-variance 
matrix)  from  the  following  data:  Observations,  standard  errors  of  the 
observations,  a  cofactor  matrix  for  each  observation  set,  and  approximate 
values  of  the  parameters.  Correlations  between  observations  belonging 
to  different  sets  are  not  considered. 


The  subroutine  COLSB  may  be  called  after  the  final  values  of  the 
parameters  are  established.  It  computes  the  corrections  (probable 
errors)  of  the  observations,  prints  and  plots  error  distributions, 
prints  identifications  of  sets  with  large  corrections  and  carries  out 
some  numerical  controls  of  the  computations. 

Sample  Problem 

The  relationship  between  x^  and  x^  may  be  linear: 

F(x^,X2  ;  ~  ^2  ”  ^2  ~  ^ 
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Frx;^)  =6? 


The  obsei-ved  points  (x,  ,x^)  may  have  different  accuracies  and 

J.  d 

there  may  be  some  correlation  between  the  coordinates  observed  so  that 
the  error  ellipses  of  the  points  observed  have  different  sizes  and 
directions,  (Such  correlations  can  be  caused  by  the  observation 
technique,  for  instance,)  COLSA  will  assign  different  weights  to  the 
points  in  accordance  with  their  error  ellipses.  The  resulting  parameters 
y^  and  y^  will  be  furnished  together  with  their  standard  errors  and 

cofactor  matrix,  (y^^  and  y^  are  correlated  because  both  are  computed 

using  the  same  data,) 

Data  for  COLSA  and  COLSB 


All  data  are  assumed  to  be  in  the  core  memory.  The  dimensions  of 
the  corresponding  arrays  are  as  follows; 

Obseivations;  X( 5,1000)  -  n  ^  5j  r  ^  1000 

Standard  errors  of  imit  weight:  ERZX(lOOO)  -  r  ^  1000 

Cofactor  matrices:  RC0R(5, 5,1000)  -  n  ^  5j  r  ^  1000 

Approximations  of  parameters:  YCAP(lO)  -  m  ^  10 

Alphanimieric  identifications  of  observed  sets,  consisting  of 
two  10 -letter  words  for  each  set:  IDER  (lOOO) 

TIFIC  (1000)  -  r  £  1000 
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In  addition  to  the  data,  a  subroutine,  say  FU,  to  evaluate  the 
function  F(x^,...,x  ;  y-,  .  •  •  •  .y  )  Is  needed.  The  routine  must  furnish 

the  value  of  F  as  well  as  the  n  +  m  values  of  the  first  partial 
derivatives  BF/dx  and  BF/By.  The  calling  of  COLSA  and  COLSB  is 
described  below. 

Calling  COLSA 

GALL  colsa(n,m,fu,ycap,ylow,ery,erzero,q,nr,x,erzx,rcor,ntot,eta,iden, 

TIFIC) 

In  the  following  list,  I  =  input  argument  provided  by  calling 
routine  and  R  =  result  set  by  COLSA. 

I-N  =  nimber  of  observations  in  each  set  (N  ^  5) 

I-M  =  number  of  parameters  (M  ^  10 ) 

I-FU  =  name  of  subroutine  for  F(x,y) 

I-YGAP  =  approximate  values  of  the  parameters.  Dimension:  YCAP(lO) 

R-YLOW  =  YCAP  +  ETA  =  improved  parameter  values.  Dimension:  YLOW(lO) 

R-ERY  =  standard  errors  of  YLOW.  Dimension:  ERY(lO) 

R-ERZERO  =  standard  error  of  weight  one  (associated  with  cofactor 
matrix  Q,) 

R-Q  =  cofactor  matrix  of  YLOW.  Dimension:  Q(l0,10) 

I-NR  =  n-umber  of  sets  observed  (NR  ^  1000) 

I-X  =  observations.  Dimension:  X(5,1000) 

I-ERZX  =  standard  error  associated  with  each  set.  Dimension:  ERZX(lOOO) 

I-RCQR  =  cofactor  matrices  of  observation  sets.  Dimension:  RCQR(5, 5^1000) 

R-NTOT  =  number  of  valid  observation  sets.  Normally  NT0T=NR;  see 
description  of  FU  -  subroutine. 

R-ETA  =  corrections  of  YCAP.  Dimension:  ETA(lO) 

I-IDEN  _  alphanumeric  identifications  of  observation  sets. 

I-TIFIC  ”  Dimension:  IDEN(lOOO),  TIFIC(lOOO) 

The  results  of  COLSA  are  printed  in  a  self -explaining  way.  The 
following  symbols  are  used  by  the  routines  for  printouts: 
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X 

X  +  KSI 

Y 

Y  +  ETA 
G 

F(X,Y  +  ETA) 


observations 

corrected  observations 

approximations  of  parameters 

corrected  parameters 

weight  of  observation  set 

constraint  function  with  corresponding 
arguments  (Equation  (l)). 


There  is  no  checking  by  COLSA  whether  N  or  M  exceeds  5  or  10, 
respectively.  Also,  the  value  of  NR  is  not  checked  to  see  if  it 
satisfies  M  <  NR  ^  1000. 

Programming  F(x,y) 

The  constraint  fimction  is  entered  by  COLSA  and  COLSB  as  follows 
CALL  FU(X,Y,FVA,A,B,L,NMES)  where 

I-X  =  observations.  Dimension:  X(5^1000) 

I-Y  =  parameter  vector.  Dimension;  Y(10) 

R-FVA  =  F(x^, . . . ;  y^, . . .y^) 

R-A  =  9F/9x.  Dimension:  A(5) 

R-B  =  3F/dy.  Dimension:  B(10) 

I-L  =  number  of  set.  The  arguments  x, ,...,x  are: 

1  n 

x^  =  X(1,L),  X2  =  X(2,L),...,x^  =  X(N,L) 

R-NMES  =  error  message  indicator.  NMES  =  1  if  the  function  or 
its  derivatives  cannot  be  computed  with  the 
arguments  given.  If  this  happens,  COLSA  will  print 
a  comment  and  the  identification  of  the  corresponding 
set  and  reduce  the  total  number  of  valid  sets, 

NTOT  by  one. 


Calling  COLSB 

CALL  COLSB(N,M,FU,YCAP,yLOW,ERY,ERZERO,Q,NR,X,ERZX,RCaR,ETA,rDEN,TIFIC) 

The  arguments  of  COLSB  are  the  same  as  those  of  COLSA,  except  NTOT 
is  not  an  argument  of  COLSB.  In  contrast  to  COLSA,  all  COLSB  arguments 
are  of  input  type.  Normally  they  will  have  the  values  furnished  by  the 
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last  call  of  COLSA.  All  output  of  COLSB  is  on  the  printer  (in 
self -explaining  manner)  and  plotter  tape.  The  plotted  information 
includes  the  frequency  distribution,  cumulative  histogram  and  probit 
diagram  of  the  weights  of  F(X,y  +  ETA).  For  comparison  the  corresponding 
number  distributions  and  normal  distributions  are  included.  The  printed 
information  includes  lists  of  sets  with  large  errors,  surveys  about  the 
corrections  and  their  distributions  and  some  numerical  tests  for 
accuracy  and  randomness. 

Subroutines  Used 


The  following  names  are  names  of  subroutines  used  by  COLSA  and 
COLSB. 

CHOLES  -  cholesky  algorithm  routine  for  solving  normal  eqiiatlons 
PLOCOB  -  plotting  routine  used  by  COLSB 

PLODIS 

FREDIS 

CUMMIS  subroutines  used  by  PLOCOB 

PRODIA 

PROBIT 

ERF(X)  -  function  routine  fiarnishing  normal  distribution 
function  (error  integral) 


PLTCCA 

PLTCCT 

PLTCCD  plotter  subroutines. 

PLTCCS 

PLTCCB 


NLPROG  (Taken  from  BRL  Memorandum  Report  No.  I958) 

In  this  section  instructions  to  use  the  NLPROG  programs  will  be 
given.  These  instructions  will  usually  be  given  as  if  POWELL  were  the 
minimizing  method.  The  necessary  modification  for  the  other  minimizing 
subroutines  will  be  noted. 

The  problem  we  wish  to  solve  is:  Minimize 

F(X)  X  =  (x^,X2,  ...,Xjj)'^ 

subject  to 

R^(X)  s  0  i=l,2,...,NIC 
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and 


R^(X)  =  0  i=NIC+l,  NIC+2, . . .,NC. 
NLEROG  act\ially  tries  to  minimize  the  function 


^MEC  1  ^NC  p 

P(X,p^)  =  F(X)  +  1/R.(X)  +  (1/p^*)  S.C) 

^NEC 

for  =  p^,  p^  =  Pj_/RATIO,  p^  =  P^/RATIO,...  until  2.i=l  ^  ®- 

NLEROG  is  really  a  set  of  subroutines  the  names  of  which  should  be 
avoided  in  other  programming.  A  list  of  the  subroutine  names  used  in 
NLEROG  with  POWEEi  follows.  The  brief  explanation  with  them  assumes 
that  X  and  the  previous  X  used  both  satisfy  all  the  inequality  constraints. 

NLEROG  -  This  subroutine  is  the  main  routine  for  this 
collection  of  subroutines. 

SUBPRO  -  A  buffer  between  NLPROG  and  the  minimizing  subroutine. 
POWELL  -  Minimizing  method. 

LINMIN  -  Univariate  minimization  scheme:  find  S  for  which 
P(XI+S*V,p)  is  a  minimum  (or  improved 
sufficiently).  This  uses  POPS  to  evaluate  P. 

POPS  -  Pinds  X  =  XI  +  S-V  and  uses  POFX  to  find  P(X,p). 

This  also  keeps  track  of  the  X  which  minimizes 
P(X^p),  say  XMIN,  and  the  functions  of  XMIN; 

R^(XMIN)  (i=l,2,. . .,NC),  P(XMIN,p),  G(XMIN,p)  and 

P(XMIN).  If  X  is  not  in  the  accepted  region 

P(X,p)  is  set  to  10^^^. 

POPX  -  Controls  computation  of  R.(X)  (i=l,2, . . . ,NC) ; 

P(X),  P(X,p),  and  ^ 

G(X,p)  =  P(X,p)  -2p  l/R^(X). 

OUTPUT  -  Prints  P(XMIN),  P(XMIN,p),  G(XMIN,p),  p,  (P-G)/2, 

XMIN,  andR.(XMIN)  ( j=l,2, . . . ,NC) . 

J 

In  addition,  the  user  must  code  a  subroutine,  EVAL(l),  which  puts  F(X) 
in  F  if  I  =  0  and  Rj(X)  in  R(i)  if  I  >  0.  If  we  use  HOOKE  instead  of 
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POVffiLL,  LINMIN  is  dropped  from  the  package.  If  FI£TCH  is  used  in  place 
of  POWELL  the  subroutine  DPBTDX  is  added  and  the  coder  must  add  the 
subroutine  DERIV(l)  where  if  I  =  0,  DRDX(j)  =  aF(X)/9Xj  ( j=l, 2, . . . ,N) 

and  if  I  >  0,  DRDX(j)  =  3R^(X)/9x.  ( j=l,2, , . , ,n) .  If  NEWTON  is  used 

in  place  of  POWELL,  the  DERIV(l)  subroutine  must  also  compute  second 

derivatives:  If  I  =  0,  DDR(j,k)  =  d^F(X)/SXj9Xj^  ( j=l, 2, .  . .  ,N; 

k=l,2, .  .  .  ,N)  .  If  I  >  0,  DDR(j,k)  =  d^j(X)/aXjaXj^  ( j=l, 2, .  .  .  ,N; 

k=l,2,...,N). 

The  original  plan  was  to  code  up  different  minimization  schemes 
and  change  only  SUBPRO  to  use  them.  Actually  the  main  control 
subroutine  NLPROG  is  slightly  different  for  each  of  the  four 
minimization  schemes  and  the  LINMIN  for  POWELL  is  different  than  that 
used  with  FIETCH  and  NEWTON. 

Blank  COMMON  is  used  to  identify  all  the  numbers  which  are  used, 
supplied,  or  needed  by  the  MAIN  program,  EVAL(l),  or  DERr/(l).  One 
labeled  COMMON,  called  SUB,  is  used  by  NLPROG.  The  blank  common  list 
for  NEWTON  is  (N,NC,NIC,X(100) ,EPS(l00),Q(l00) ,RH0, RATIO, THETA, NREPET, 
NRHORP , F , R ( 200 ) , IPR OP , DRDX ( 100 ) , DDR ( 100 , 100 ) ) . 

N  -  number  of  variables.  N  ^  100. 

NC  -  number  of  constraints.  0  s  nC  ^  200. 

NIC  -  number  of  inequality  constraints.  NIC  £  NC  s  200. 

X(lOO)  -  variables.  An  initial  guess  must  be  put  here.  The 

current  value  of  X  is  here  at  other  times. 

EPS(lOO)  -  acceptable  absolute  error  in  variables.  NLPROG 
may  not  attain  this  accuracy.  It  does  check 
that  a  change  of  EPS(i)  in  any  one  x(i)  will  not 
reduce  the  final  P(X,p).  EPS(i)  >  0. 

Q(100)  -  initial  search  steps  (not  too  critical) .  The 

approximate  errors  in  the  initial  x.  are 
recommended. 

RHO  -  p.  control  parameter,  p  >  O. 

RATIO  -  =  p^/RATIO.  RATIO  >  1  (Recommend 

4  ^  RATIO  5  16.) 

THETA  -  9.  final  convergence  parameter. 

NREPET  -  number  of  times  to  repeat  minimization.  (Optional) 
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NRHORP  -  number  of  p's  to  repeat  for  each  repetition  of 
minimization.  If  is  the  value  of  p  when 

P  l/Ri(x)  <  0  the  minimization  is  restarted 

with  p  =  Pt  .  The  purpose  of 

Lf 

this  is  to  push  the  X  which  satisfies  P(X,p)  away 
from  the  boundary  of  the  region  {x|r^(X)  >  0, 

i=l,2, . . . ,NIC} .  However  if  there  are  no 
inequality  constraints,  only  equality  constraints, 
a  negative  or  zero  NRHORP  would  reduce  p  and  force 
closer  satisfaction  of  the  equality  constraints . 

F  -  function  to  be  minimized. 

R(200)  -  constraints.  (The  NIC  inequality  constraints  must 

be  stored  first.) 

IPROP  -  controls  use  of  OUTPUT. 

0  -  no  prints  (except  errors) 

1  -  print  initial  values  and  the  solution  of  each 

subproblem. 

2  -  Print  each  cycle  result 

5  -  Print  each  linear  minimum.  (iPROP  =  1  is 

recommended.  This  gives  a  fair  outline  of 
the  course  of  the  problem. ) 


DRDX(lOO)-  Partial  derivatives  of  F  or  R^  with  respect  to  x^ 
(This  array  is  not  needed  for  POWELL  or  HOOKE.) 

DDR(100,100) -  Second  partial  derivatives  of  F  or  R^  with  respect 

to  x^  and  x^.  This  matrix  is  symmetric  but  must 

be  completely  filled  by  DERIV(l)  (This  array  is 
needed  only  if  NEWTON  is  used. ) 


The  quantities  starting  with  N  or  I  are  integers.  The  rest  are 
floating  point  numbers.  This  set  of  blank  COMMON  must  precede  the  users 
programs . 

Before  calling  NLPROG,  the  user  must  store  N,  NC,  NIC,  and  initial 
values  of  X(l),  W(l),  and  EPS(l)  (l=l,2, . . . ,n) .  He  must  also  store 
RHO,  RATIO,  THETA,  IIROP,  and  if  desired,  N^PET  and  NRHORP.  He  must 
supply  the  subroutine  EVAL(I)  and,  if  FLETCH  or  NEWTON  are  used,  the 
user  must  also  supply  the  appropriate  DEEIV(l)  subroutine. 
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When  NLH?OG  returns  control  to  the  main  progreim,  X(l)  1=1. 2,..., N 
contains  the  solution  RHO  has  the  final  value  of  p,  F  =  F(x),  and 
R(J)  =  R,(X)  J=l,2, . . . ,NC.  The  numbers  in  N,NC, NIC, EPS(l) , RATIO, THETA, 

NRHORP,  and  IPROP  are  not  changed  by  NLPROG. 


Comments 

The  need  for  some  preliminary  analysis  of  minimization  problems 
cannot  be  overstressed.  All  too  frequently  we  mechanically  prepare  a 
program  and  discover,  while  trying  to  code  check  the  program,  that  there 
is  an  obvious  solution;  or  worse  that  there  is  no  finite  solution. 
Further,  the  preliminary  analysis  should  help  select  reasonable  starting 
values.  If  the  initial  X  does  not  satisfy  the  inequality  constraints 
NLIROG  will  search  for  an  X  that  does.  However,  this  X  may  be  so  remote 
from  the  desired  solution  that  the  program  will  take  a  long  time  to 
reach  the  solution. 
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